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THE AUTOMATIC DESIGN OF STRUCTURAL FRAMES 


By R. K. LIVESLEY 


(The Engineering Laboratories, University of Cambridge) 
[Received 30 June 1955] 


SUMMARY 


This paper considers the problem of finding the lightest structural frame of given 
geometrical form which will support a given set of loads. Following the development 
of a geometrical analogue and an iterative method of solution, an analytic technique 
is presented which gives an exact solution in any particular numerical case. The 
method is very suitable for use on an electronic computer, and a brief description 
of programmes developed for the Manchester University machine is included. It is 
shown that with slight modifications these programmes can also be used to determine 
the collapse loading of a given frame. 


1. Introduction 

THE problem considered in this paper may be stated as follows: Given a 
set of static concentrated loads, acting at certain fixed points of a rigid- 
jointed plane frame of prescribed geometrical form, how should the cross- 
sectional dimensions of the members be chosen to produce the lightest 
possible frame capable of carrying the loads? The members are required 
to be straight and of constant cross-section throughout their length, and 
the term ‘frame’ is used to denote a structure which resists deformation 
entirely by bending action within its members. 

To determine whether a frame will support the loads applied to it the 
theory of plastic collapse will be used (1). This theory applies to structures 
made of a ductile material, and assumes that if the curvature of a member 
becomes infinitely large the bending moment tends to a maximum value, 
called the fully plastic moment, which depends only on the section dimen- 
sions. If this moment is attained at a particular cross-section, a plastic 
hinge is said to have formed, and at such a hinge a finite change of slope 
can occur in the member. The magnitude of this change is independent 
of the moment, while the direction of the moment is always such as to 
oppose further rotation. Provided that instability and shear effects can 
be neglected, a structure will be on the point of collapse if sufficient plastic 
hinges exist to transform it into a mechanism, since under these conditions 
deformations can increase indefinitely without further increase of load. 


This criterion of failure implies that in any given problem a ‘design’ can 


be specified by giving particular values to the fully plastic moments of 
all the members of a frame. To simplify the analysis it will be assumed 


{Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 3 (1956)] 
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258 R. K. LIVESLEY 


that the fully plastic moment .@ of a member is proportional to its cross. 


sectional area.t This implies that the total weight of a structure is a linear 
function of the fully plastic moments, and the problem is to minimize this 
function subject to certain constraining conditions. These conditions are 
imposed by the fact that the moments in the structure must be in equili- 
brium with the known external loads, and must be less than or equal to 
the fully plastic moments of the members in which they occur. Mathe- 
matically, therefore, the problem is one of minimizing a linear function of 
several variables subject to a set of linear inequalities. This is the basic 
problem of linear programming, so that any general method of designing 
minimum weight structures may be of interest in other contexts. 

Methods of minimum weight design under the conditions outlined above 
have been discussed by Foulkes (2, 3) and Heyman (5). The present paper 
develops an alternative system of analysis which provides an exact solution 
of the general problem, and which is particularly suitable for use on an 
electronic digital computer. The paper describes briefly two general pro- 
grammes which have been constructed for the machine at Manchester 
University, and discusses some of the computational problems which arose 
in their development. No detailed account of automatic computing is 
included, but examples are given to demonstrate the speed with which the 
machine can solve particular design problems. 


2. Notation and mathematical formulation 

If a structure is acted on by concentrated forces only, the maximum 
bending moment in each member occurs either at an end or at a point 
where a load is applied. In a given problem, therefore, there is a fixed set 
of points in the structure at which plastic hinges can occur, and the moments 
at these points will be denoted by M,, where in general « runs from | to n. 
This set of moments can be divided into smaller sets or ‘groups’, each group 
corresponding to points at which, by the formulation of the problem, the 
member cross-sections are known to be equal. In most cases a group merely 
comprises all the moments associated with a particular member, but in 
certain problems two or more members might be required to have equal 
cross-sections, and the moments in all these would then be included in a 
single group. Since the geometrical form of the structure is assumed to be 
fixed, there will be a certain known length associated with each member 
or group. It is convenient to include these lengths in the analysis by 
introducing weighted moments X,, derived by multiplying each moment 
M, by the length / associated with the group to which that moment belongs. 

Ifa frame has r redundancies, the equations of equilibrium may be written 


+ The accuracy of this assumption is discussed elsewhere (3, 4). 
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in a form which expresses the moments , in terms of r arbitrarily-chosen 
moments m. (7 a r), 


M a,,m,+MO. (1) 


x x 
Multiplying (1) by the appropriate lengths gives a similar expression for 
the weighted moments X , 


X A ,;mt+X. (2) 


x Xt 


The matrices a,;, A,; are n <r rectangular matrices depending only on the 
geometry of the frame, while the columns M®, X depend also on the 
applied loading. For the beam shown in Fig. 1, for instance, plastic hinges 
can form only at the points 1, 2,..., 5. If the moments at the points of 


support B and C are chosen as redundant moments, elementary statics 


gives 
Group | ([X, l 0] [m,]+16 
] 2 \| X, 2 0 a lo 
Group 2 | As : . i (3) 
l 4 | < ~ ye . 
X 0 4 0 


Any set of values VM, or X, which satisfies equations (1) or (2) is termed 


a statically admissible set. 
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For any given moment system it is possible to associate with each value 
of « a parameter 4, defined as the suffix of the moment of maximum 
modulus in the group to which M, belongs. Since the weighted moments 
X, within a particular group are constant multiples of the corresponding 
moments V,, this definition implies that |X;| > |X,| for all values of «. 
The set of numbers & form a sub-set of the m suffixes a, each value of @ 
being associated with a different group. In the analysis which follows, the 
symbols « and 7 will be used as general suffixes, repetition implying sum- 
mation, while other letter suffixes will denote particular variables or 


equations. The symbo will be used to denote s ation over the 
juat The symbol 5 will I 1 to denote summat t! 


sub-set 4, the components of this set being defined as described above. 
The assumption that the fully plastic moment .@ of a member is pro- 


portional to its cross-sectional area implies that the weight of a structure 
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is proportional to the quantity }.@/, summed for all the individual 
members or groups. For any arbitrary set of-moments m,, equation (1) 
defines a statically admissible set M,, from which the values of « forming 
the a-set can be picked out by inspection. These values specify the point 
in each group at which the moment is greatest, and hence the points in 
the structure where the particular moment system M, is most likely to 
induce plastic hinges. The fully plastic moments of any structural design 
capable of supporting the moments M, must clearly satisfy inequalities 
of the type > |M;|, one such inequality being associated with each 
group. The lightest of these structures will be the one in which equality 
holds in every case, i.e. 

A= |M;|, or al = |X;|, (4) 
and this structure will be termed a critical design. It will have at least one 
plastic hinge in each group of points, and the weight will be proportional 
to the function G, where, from (4), 

G > X =|. (5) 
Substituting for X; from (2) gives 
G = > |Ax;m;4 x 
1 


xi i X 


= ¥ (A;;m,;+X2)sgn X;. (6) 


The equalities represented by equations (4) associate with each set of 
redundants m, a critical design which can just sustain the moment system 
M, appropriate to that set. It follows from the lower bound theorem of 
plastic collapse that all such designs will in fact sustain the applied loads.t+ 
Since by a suitable choice of values of the redundant moments it is possible 
to generate any statically admissible moment system, it follows that if a 
set of moments m; can be found which minimizes the function G(m;) defined 
by (6), then the associated critical design given by (4) will be the least 
weight structure required. 


3. A geometrical analogue 

Although not strictly necessary to the analysis, it is useful to adopt the 
language and concepts of, geometry in discussing the function G(m;,). The 
geometrical representation developed here is due to Jennings (4). It uses 
the redundant moments m; as coordinate variables and therefore differs 
from the analogue suggested by Foulkes (2). 

If the redundants m, are considered as defining an r-dimensional Euclidean 

+ This theorem may be stated as follows: ‘If for a given structure acted on by given loads 
it is possible to find a statically admissible moment system such that the fully plastu 


moment of no member is exceeded, then the loading is not greater than the collapse load 
of the structure.’ See, for instance, Greenberg and Prager (6). 
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lual | space, equation (1) associates with each point in that space a moment 
(1) | system ,, and hence a particular critical design and value of G. The 
ning | complete space may be divided into ‘regions’, a region being defined as 
oint 
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the set of all points for which the moments of maximum modulus M, occur 
at certain fixed places in the structure and have constant sign. Throughout 
a region the sub-set & remains unchanged, so that the function G defined 
by (6) is linear in the variables m;. In general, the regions are convex 
polyhedra, their boundaries being hyper-planes on which either two 
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moments within a group have equal modulus (greater in value than any 
other moment in the group) or some moment of maximum modulus changes 
sign. Movement across a hyper-plane from one region to another implies 
a change in either the sign or the position of a moment of maximum 
modulus, and hence in the coefficients of the linear function G(m;). In 
the example shown in Fig. 1, for instance, the m;-space is two-dimensional, 
and is divided into regions as shown in Fig. 2. Lines of constant G are 
shown on the figure, together with the positions and signs of the dominant 
moments in the various regions. 

Before discussing methods of finding the minimum of the weight function 
G, it is necessary to prove that this function has, in fact, a unique minimum 
value. Since, from physical considerations @ is always positive, and is 
continuous throughout the m,-space, it is sufficient to prove that, for any 
constant A for which the surface G = ) exists, the surface is bounded and 
convex, and that G > A at any point external to it. 

From equation (5) it follows that ¥ |X; A on the surface, and hence 


; 
that X;| <A. Since the weighted moments X, are linear functions of the 
variables m,;, the latter must also be finite, and the surface is therefore 
bounded. 
To prove that the surface is convex, consider the hyperplane 
A= > (Agim,+X)sgn XF, (7) 
rm 

where jf indicates the set 4 associated with some region R which the surface 
intersects, and X* indicates the set of quantities X;, at some point in R. 
Let 7 indicate the components of the d-set at some point S on the hyper- 
plane, where S lies in a region R’ different from R. Then in R the hyper- 
plane defined by (7) is coincident with the surface G = A, while at S, 
X;| > |X,|, with inequality in at least one case. It follows from (5) that 


G, = > |X;| > A. Thus G > X at all points on each hyper-plane forming 


5 
part of the surface G = 4, so that all the conditions required for uniqueness 
are fulfilled. 

From the linear nature of the function G, it is apparent that the minimum 
will always be attained at the vertex of some region, where r hyper-planes 
intersect. If two or more vertices give rise to critical designs of equal 
minimum weight, then all points belonging to the convex set generated 
by these vertices will correspond to designs of the same weight, but other- 
wise the minimum will occur at a unique point. In the example shown in 
Fig. 2, the minimum vertex is the point #, corresponding to the solutior 
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4, An iterative method of steepest descents 


In the geometrical analogue described above, each vertex is defined by 
moment-equalities, each equality relating two moments belonging to the 
same group. The problem is therefore to find the particular set of equalities 
which defines the vertex at which G(m,) is least, and the set of weighted 
moments X , appropriate to that vertex. 

Both the methods described in this paper attain the minimum vertex 
by tracing a path through m,-space from some initial trial point m$”. In 
each case the path consists of a series of discrete steps, and it is convenient, 
before each new step is taken, to make the trial point which has just been 
reached into the origin of coordinates. For any trial point m equation 
2) defines a ‘trial solution’ X, so that if new independent variables 2; are 


defined with the point m\? as origin, equation (2) can be written 


F ° ™(k) » 
X, = A,,2%;,+XY. (8) 
In the same way equation (6) becomes 
’ r(k , 
G = > (A;;a,+Xz’)sgn X;. (9) 


In the analysis which follows, the variable x; will be used in the sense 
indicated above, representing a displacement from the current trial point. 
Since the solution of the problem can be obtained direct from the last term 
of the sequence X®, X,.... X..., it is unnecessary to keep a record of 
the absolute coordinates m; of each point. When a new solution X%*" is 
found, the shift of origin is achieved merely by substituting it for the old 
solution X“ in (8). 


The simplest method of tracing a path towards the minimum vertex is 


by moving a suitable distance in the ‘direction of steepest descent’ at each 


trial point. Differentiating (9) at the point m\” gives 
oG > (Kk) 
> Az,;sgn Xz” = p;, say, (10) 
CA 


where the set & is determined by inspection of the vector X”. It is clear 


that the components of the gradient vector p; are constant throughout the 


region containing m\). If @ is some positive number, then movement to a 
new point m“*+» defined by 
v m'*+1)_m\*) — —Op,, 


represents a ‘step’ @p,; in the direction of maximum decrease of G. The 
trial solution associated with the new point is given by 

"(k+1 ”(k) 

xt XX —OA |; Dis (11) 
and the process can clearly be repeated indefinitely. The path will eventually 
reach a state of random variation about the minimum vertex, the speed of 
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convergence and the accuracy of the solution depending on the parameter 6 
which controls the step length. 

The method outlined above is extremely simple to programme for an 
electronic computer, and has been extensively tested on the Manchester 
machine by Jennings (4). In the simplest of these experimental programmes 
the sequence of operations is as follows: 


(1) The computer is supplied with the matrix A,;, the vector X, and 
the coordinates of the first trial point m$”. In addition, a number is 
fed into the machine whose binary digits specify the division of the 
moments into groups. 

(2) The trial solution X{ is formed according to equation (2), and this 
is substituted for the original vector X, thus shifting the origin to 
the point m\”. 

(3) From the trial solution the components of the sub-set & are deter- 
mined, and the vector p; calculated according to (10). 

(4) The new trial solution is calculated from (11), the parameter @ being 
chosen by the operator and fed to the computer by means of switches 
on the machine’s control panel. 

(5) The new solution and the value of G are printed if required, and the 
programme returns to step (3), with X{ replaced by X®. After k 
cycles the solution X‘ is thus replaced by X%+, 


This iterative scheme proceeds automatically, apart from the manual 
setting of 0, which may be altered at any time during the calculation. The 
time required for a single step varies with r and n—in the example shown 
in Fig. 2 each step took about 2 seconds. 

The chief difficulty experienced with the programme was in adjusting 
the parameter @, and in deciding when a solution reasonably close to the 
minimum had been reached. It is clearly desirable to commence with a 
fairly large step and to decrease 6 as the minimum is approached, but it 
was found in practice that a premature reduction often resulted in a false 
limit point being attained. An automatic method of controlling the step 
size would have been preferable, but no criterion was found which could 
be applied in all cases. , 

Another disadvantage was the time factor. Experience showed that in 
many problems the steps were so misdirected that many hundreds were 
required to reach a satisfactory solution. It is apparent that when the 
weight function G has gradient vectors in two adjacent regions which are 
almost in opposite directions, the trial point may alternate continuously 
between the two regions without making much progress in the right 
direction. This situation may then deceive the machine operator into 
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ré@ \ thinking that the step length should be reduced, in which case a false 
| solution will be reached. 
an | These factors render the iterative approach unsatisfactory except in 
ter simple problems with only two or three redundancies. An alternative 
1es method developed by the present author which avoids these convergence 
difficulties will now be described. 

nd 

‘ 5. A modified method of steepest descents 

he As before, the method requires a trial solution X) derived from an initial 


trial point m\”. Subsequent work falls into two distinct stages. If R is the 
region to which the initial point belongs, then the first part of the analysis 
t. » is concerned with finding the solution associated with some vertex of R. 
This corresponds to the ‘basic feasible solution’ of linear programming. 
1 The second part of the work involves movement from vertex to vertex 
until the minimum is attained—a procedure very similar to that which 
1g occurs in the ‘simplex’ method (7). 


5.1. Attaining a vertex 

In the iterative scheme described in § 4 the path always follows the line 
of steepest descent appropriate to the region containing the current trial 
point. In the method now to be described, the path follows the line of 


: steepest descent until it reaches a bounding hyper-plane of the initial 
al region R. It is then constrained to lie within this boundary, following a 
le modified line of steepest descent until it intersects another hyper-plane, 
n> after which it moves through the sub-space common to both. Each inter- 


section imposes one more constraint on further movement, a vertex of R 
being reached when a total of r constraints have been imposed. 


e At the initial trial point m\? equation (10) defines the gradient vector p; 
a associated with the region R. The line of steepest descent is given by 
it 
? . 9 

‘ x; = —Op,, (12) 
p where @ is now considered as a variable which may take any positive value. 
d | Writing P, A; p;, Substituting for x; from (12), and taking the point 

m;” as origin in equation (8) gives 
n ; 

4 D 7 (1) y. 

e 5 xX, OP. + Ay”. (13) 
p This equation defines the varietion in each X, with movement down the 
p line of steepest descent. As 6 is increased from zero, the point x; defined 
a by (12) traces a straight line through R and eventually this line must 
intersect a hyper-plane separating R from an adjacent region R’. As 


» | mentioned in § 3, a hyper-plane is defined by an equality between two 
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moments in a group. Let the one dividing R from R’ be defined by the 
equality Y ei 


,=4+X,, (14) 
where p and v are particular values of «a, and where 


X X;| > |X in R, 


/ F . 
A, X,;| > |X,| in R. 

In other words, the quantity X,, is the dominant one of its group in R, 

while in R’ it is deposed by X,. The equality (14) will of course have a 

definite positive or negative sign in any particular problem. 

The complete line defined by (12) intersects a large number of hyper- 
planes defined by equations similar to (14). To find the correct equality 
it is necessary to determine, for the set & appropriate to R, the value of 6 
associated with each of the equalities 


X,;=1+tX (a ~ a), 


: EA, (15) 
and then to select the values of « and &, for which @ takes its smallest 
positive value. From (13) and (15) 
(x 1 xX) 
xX 
| P, ¥ i. 


where the sign depends on the sign of equality (15). The possibility of a 


@ (a ~ a), (16) 


positive or a negative equality must, in general, be considered for each 
value of a, but it will be shown later that a simple test can be devised to 
determine which sign is important. For the moment it will be assumed 
that the smallest positive value of @ given by (16) is 6,, where yu, v and the 
sign of the equality in (14) are known. 

The next stage of the analysis involves a linear transformation of 
equation (8). The first part of this transformation consists of a shift of 
origin to the point 2; 6, p; (i.e. the point at which the line of steepest 
descent intersects the boundary of R given by (14)). From (13) the trial 
solution X® associated with this point is given by 

A® = 0, P.+-X, (17) 
so that referring the variables x; to the new point, equation (8) becomes 
X, = Ay t%+X. (18) 
The second part of the transformation consists of a change of variable. 


A new variable y, is defined by the equation 


Yr : Ai I A “Ne 0, 
ol Fe 


pl 


(19) 


where / is the value of 7 for which the term (A,,;- A,;) has largest modulus, 
and the sign depends on the sign of equality (14). This method of choosing 
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the suffix / is not strictly necessary to the analysis, but is a convenient 
way of keeping the equations well-conditioned. It corresponds to the 
process of taking the largest pivot in the elimination method of solving 
simultaneous equations. Multiplying (19) by A,, and adding to (18) gives 


? A, it A,j;' . 7 (2) 2 
X A, | i "' A uf +A YytX. (20) 
fe lV Seog 

It will be seen that the coefficients of 2, in (20) are identically zero, so that 
the operation has the effect of replacing the variable x, by the new variable 
y,. Equation (20) may be written 


X, = Ay 2;+X, (21) 
where 2x, indicates the set of variables 2,, %9,..., Y,-.-, 2,, and A}, indicates 
the matrix obtained by subtracting the appropriate multiples of A,, from 
the other columns of the original matrix A ,;, asin equation (20). The column 
A,, itself remains unaltered, although it is now associated with the new 
variable y 


It is easy to verify from (20) that 


A'.=+A', (64D), (22) 


where the sign is the same as the sign of the equality (14), and it follows 
that 


X, FX, = (AwFAny (23) 
all terms in the other variables dropping out. Thus any variation from 
the new trial point in which y, is kept equal to zero does not affect the 
equality (14), and therefore corresponds to a movement in the hyper- 
plane defined by that equation. In particular a new direction of steepest 


descent p; defined according to the equations 


OG “ r(92 . 

pi = =TAgsgenX” (AD, ) 
ax, = (24) 
0 (2 »,J 


lies in the hyper-plane, where the same set & may be used as before, since 
the point just reached can still be regarded as belonging to R. From the 
new gradient components a new set of parameters P’, = — Aj, p; may be 
found and the whole process repeated. In this way a new trial solution 
X is reached and a new variable y,, introduced. The choice of y,, is made 
in the same way as before, with the restriction that m must be different 
from /. Any variation in which both y, and y,, are kept equal to zero will 
not disturb either of the equalities which have been associated with these 
variables, and therefore corresponds to movement within the sub-space 
common to two of the hyper-planes bounding R. 


Each repetition of the analysis results in one more equality being set 
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up among the weighted moments X,, and one more restriction being 
imposed on variation in the variables x;. As each new variable y,, y,,,... 
is introduced, so the number of non-zero gradient components in equation 
(24) is reduced, until finally only one of the original variables remains, and 
(r—1) equalities similar to (14) have been set up. This implies that a line 
joining two vertices has been reached, and the next application of the 
process produces a movement along the line to one of these vertices. The 
r equalities similar to (14) now determine the coordinates of the vertex 
completely, and in fact the solution has already been obtained in the form 
of the vector X*" generated by repeated application of (17). 


5.2. Finding the minimum vertex 

The series of transformations described in § 5.1 achieves two results. In 
the first place it alters the initial trial solution X® to a solution X%+)) 
associated with some vertex of R. The components of the vector X‘+!) 
satisfy r equalities of the type |X | = |X ,| (a # «), where & indicates the 
set appropriate to R. In the example shown in Fig. 2, for instance, two 
cycles of the analysis result in the vertex C being reached from the initial 
trial point A, via the path ABC shown. In the second place the original 
matrix A,; appearing in (2) is transformed into a new matrix A/‘,, asso- 
ciated with a new set of variables y;, so that each such variable is associated 
with relaxation of only one of the equalities which together define the 
vertex. 

If the vertex which has been reached is the one for which G attains its 
minimum, then any further movement will increase G. If the minimum 
has not been attained, however, there must be at least one line connecting 
the vertex to another of less weight. If such a line can be found, the analysis 
already developed can be used to determine a new vertex. If several lines 
exist, the best results will be achieved by choosing the line on which G 
decreases most rapidly. This is not essential, however as the process will 
always terminate if any of the lines on which G decreases is chosen at each 
Stage. 

In attaining the first vertex the path lies either in R or on its boundaries. 
The set & therefore remains the same as it was at the original trial point 
m\». In considering further movement it is convenient to distinguish 
between variations which lie along the boundaries of R and those which 
do not. 

5.2.1. Variation along a boundary of R. The original equation (8) has 
now been transformed into 


Xx _ Ayiyit Xf ; *, (25) 
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Consider the variable y, associated with the equality X,, = +X,, where 
X,, is the dominant quantity of its group in R. If only y, is varied, (25) be- 


comes , ’ 1 Wire 92 
X AyyytXer, (26) 


x 
and from previous analysis such a variation will not upset the (r—1) 
equalities associated with the other variables. These (r—1) equalities 
determine a line passing through the vertex under consideration, so that 
variation in y, represents a movement along this line, and a relaxation 
of the related equality (14). Differentiating (26) at the vertex gives 


c 


to 
~! 


(|X,,|—|X,|) = Aj,sgn X0+Y —A}, sgn XPt”. (: 
< Yi , 


Writing B,,, for the right-hand side of equation (27), it is clear that the 
variation in y, for which |X,,| > |X,| is that which has the same sign as 
B,,. This variation implies no change in the a-set, and represents move- 
ment along the line in the direction in which it forms a boundary of R. 

[t is apparent that the column of coefficients A), in (26) plays exactly 
the same role as the column P, in equation (13). Equation (26) gives 


0G 


‘ Pee y(r+1) 
> A y Sgn X * ’ 
CY) Xx 
c G > = , ' r(r+1) ’ 9 
and hence sen B., > Ax sgn XZ sgn B.,,, (28) 
CY, ¥ 


where, as before, 4 represents the set appropriate to R. Equation (28) gives 
the differential coefficient of G with respect to a variation y,sgn B,,,, Le. 
the direction for which the line in question forms a bounding line of R. 
If this expression is negative, then a movement from the vertex to a trial 


int where 
- om Y; @ésgn By (6 > 0) 


and the other variables remain zero will result in a decrease in G. Multiply- 
ing the column A‘, by sgn B,,, results in a situation in which positive 
variation in y, will reduce G. 

5.2.2. Variation away from R. If the minimum vertex is not a vertex 
of R, then at some stage a point will be reached where a decrease in G can 
only be achieved by a movement along a line leading away from R. Such 
a movement involves a change in the position of a moment of maximum 
modulus and hence in the set 4. Two cases must be distinguished. 

(a) No other equality within the group. If the equality X, = +, asso- 
ciated with the variable y, is the only one in its group, then a variation in 
which X, becomes the dominant quantity will not affect any other equali- 


ties, since these do not contain X,. In such circumstances a variation in 


a 


which |X, X,,| represents movement along a line in the opposite 


jt 
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direction to that discussed in § 5.2.1, i.e. a movement in which y, has the 
same sign as —B,,. The differential coefficient of G in this direction is 
oG 
— — sgn B.,,, (29) 
oN 
where the summation giving @G/éy, now includes the term A‘, sgn X("+), 
associated with the new dominant quantity X,, in place of the term 
A),sgn X\7*). Expression (29) may be written 


0G 
(ar) ; B,.| en Bays (30) 
where (€G/ey,), indicates the summation previously carried out by using 
the a-set appropriate to R. If expression (30) is negative, then, as before, 
multiplication of the column Aj, by —sgn B,,, produces a situation in 
which positive variation of y, reduces G. 

(b) Several equalities in one group. If there is only one equality in a 
group, then both positive and negative variations in the associated variable 
are permissible. In geometrical terms, the line associated with the variation 
extends on both sides of the vertex, as at point C in Fig. 2. When there is 
more than one equality in a single group, however, each line will end at the 
vertex instead of continuing through it. The reason may be seen from a 
consideration of the nature of the equalities. Let there be s equalities 


x Xj, (X.j=—(X |, |X.) = (xz 


ph v 


associated with the variables y, y,,, 5. where X,, is the dominant 


Ku 
weighted moment in R for a certain group. It is possible to vary y, in 
such a way that |X| > |X,| without affecting the other equalities, but 
movement in the opposite direction renders these equalities meaningless, 
since they do not now contain the new dominant quantity X,. This 


situation occurs at points D and E£ in Fig. 2. The equations of the s lines 
lying on the boundary of R are 


X, X,| = |X) X,| = 
. it c Mi, 
X) S Xi —" x, = (31) 
Xx p X fe eee X A = 
These equations imply an additional line lying away from R 
X,| < |X,| = |X) X,| = +. (32) 


The lines defined by (31) have already been considered in § 5.2.1, while the 
implied line (32) can be dealt with as follows. 
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THE 


Consider first the equalities |X| = |X,|, |X,,| = |X)|, associated with 
the variables y,, y,,. If a new variable y,, is defined by 
Ay, +A; 
Ant Au as 
Yn Ym ' = ry Yp (33) 
“2m T 4*ym 
where the sign depends on the sign of the implied equality X, = +X), 


substituting for y,, in (25) gives 


Ay +A, ’ il : ' 
X x 1 x1 Yy A 7 ri )4 on T A xm Ym T fA vr Yr T ae ») 
“2m T “*ym 
(34) 


It can easily be verified that if A‘, represents the new column of coefficients 
of y, appearing in (34) then 

A‘, + A}, (35) 
according to the sign of the implied equality. Thus if y, is varied keeping 
all other variables including y;,, constant, the equality |X,| = |X| is not 
affected. 

In exactly the same way, by subtracting a suitable multiple of the 
column A‘, from the column A‘, a new column A“) can be derived in which 
X,| = |X,|. The 
process may be continued until all the equalities occurring in the group 


a variation of y, does not upset the equalities |X, 
have been dealt with, and the column finally obtained may be treated for 
variation along the line given by (32) in a similar manner to that described 
under (a) above. 

5.2.3. Determining the new vertex. Once a line has been found along 
which positive variation of the appropriate variable produces a decrease 
in G, the analysis developed in § 5.1 can be applied. Apart from a multi- 
plying factor the column A‘, (or the combined column in case (b) above) 
associated with y, is equivalent to the column P, associated with the 
variable @ in equation (13). Provided that the solution @ = 0 (associated 
with the vertex already found) is discarded, the smallest positive value of 
# given by (16) defines a new equality, which can be associated with the 
variable y, in the same way as before. Thus a new vertex is reached and 
the matrix A), modified by the column A‘, in the usual way. When a 
vertex is reached at which G does not decrease on any line, that vertex must 


be the required solution, and (4) then defines the associated critical design. 


6. A fully automatic programme 

Since the steps comprising the analysis presented in § 5 are all simple 
operations of linear algebra, the development of a computer programme 
to carry out these steps automatically presents few intrinsic difficulties. 


The programme developed by the author for the Manchester Computer has 
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a capacity of 32 moment points and 16 redundants, but there is no reason 
why this could not be extended. 

As in the simple iterative scheme described previously, the programme 
commences by reading the matrix A,;, the vector X, and the coordinates 
of the initial trial point m‘\”. The computer then evaluates and prints the 
initial trial solution X{”, and records the a-set appropriate to the initial 
region R. Subsequent operations follow exactly the analysis set out in § 5, 
the machine finally printing the solution associated with the minimum 
vertex, together with the value of G. The solutions associated with inter- 
mediate vertices may also be printed if required. 

The programme was conceived as a ‘standard’ programme, which could 
be applied to any frame within its capacity by practising structural 
designers. In a recent paper (Livesley and Charlton, 8) the authors 
described a similar programme for elastic- structural analysis, and dis 
cussed the general principles underlying the use of electronic computers 
in solving routine engineering problems. Since a standard programme 
may be used by many different machine operators, it is essential that it 
should work correctly when applied to any conceivable problem within its 
range. This means that it must be constructed and tested extremely care- 
fully, and that all ‘special cases’, however unlikely, must be allowed for. 
In the programme at present under discussion, a large part of the develop- 
ment time was in fact spent on making the programme function correctly 
in all situations. Much of this work was merely a matter of organizing the 
rather complicated record-keeping required, but certain points seem of 
sufficient computational interest to be worth recording here. 

The chief difficulty, which appears to occur quite generally in linear 
programming work, is due to the fact that the course of the calculation 
depends at certain stages on the equality or disparity of certain numbers. 
The section of the analysis devoted to determining the next equality 
encountered is particularly sensitive to small differences in numbers which 
are supposed to be exactly equal, and a malfunctioning of this part of the 
programme which results in the same equality being selected twice has a 


disastrous effect on the subsequent calculations. If the equality | X,,| = |X, 
is the first to be selected, for instance, then from equation (17), 
Xi | = |X? |. (36) 
On the second cycle of analysis, when the vector P’, = — A‘, p; is calculated, 
equations (22) and (24) imply 
P,,| = |P%|, (37) 


so that equation (16) gives an indeterminate value for 4 when the equality 


(14) is considered on all subsequent occasions. It is straightforward to 
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arrange that the programme ignores this value, but if rounding errors 
cause equalities (36) and (37) to be only approximate, the programme may 
mistake the equality |X, X,| for one not already attained. For this 
reason the routines are arranged so that equalities of the type represented 
by (22), (35), and (36) are always digitally exact. In the case of the matrix 
alteration specified by (20), this exact equality is achieved by calculating 
the columns of the new matrix according to the equations 


pl T A,,)A xi (A,;FA,,)A xl 
A.) t A, 


instead of the algebraically equivalent expression appearing in (20). A 
similar procedure is used in calculating the column A‘, in (34). 

As has already been mentioned, the occurrence of both positive and 
negative equalities similar to (15) must be considered when determining 
the next equality encountered. For machine: which have a relatively slow 
division time a test is useful to find which equa.ity is relevant. It is straight- 
forward to show that the sign of the expression 

P, X3)—P, X, (38) 
is the same as the sign of the equality (15) which gives the smaller value 
of 6. Thus if (38) is positive, the relevant value of 6 will be that asso- 


ciated with the equality X; = X,, ie. 
7(1) 7(1) 
A oo (X5 XY’ | 
x ) ) P 
| PP, | 
while if it is negative the only possible equality is the negative one 
Xx = 


In some cases, of course, neither equality will give a positive value of @. 
The section of the programme which tests the lines meeting at a vertex 
is arranged so that it follows the first line it comes to on which G@ has 
negative gradient. The variables y,; are scanned in turn and the appropriate 
columns of the matrix A‘, tested as described in § 5.2. It would doubtless 
be better from the point of view of convergence to test all the lines and 
choose the one for which the gradient has largest negative value, but 
successful functioning of the programme on a number of test problems 


indicates that this is not absolutely necessary. 


7. Examples of problems solved by the programme 
The following examples were among those used during the development 
and testing of the computer programmes described in this paper. They 
5092 ¥ 
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Example 1 


100 


Example 2 


function is depicted in Fig. 2. 
£ 


Input of problem tape 
Calculation time 
Printing of solution 


be the same for all problems. 


i<— /0 
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100 





<—— /0 ——> 


shown in Fig. 2. The solution time was as follows: 


50 





were originally chosen to assist in the eradication of various logical faults 
and programming errors, but they give some idea of the possibilities of the 
machine in structural design work. 


The first example is the continuous beam shown in Fig. 1, whose weight 
As a typical case, starting from the point A 
(mp = 2:75, mg = —4-25) the computer traced the path ABCDE, as 


These figures do not include the time required to feed the programme 
tape into the machine. This extra time (about 45 seconcs) will, of course, 
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equations of equilibrium are as follows: 
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This structure is shown in Fig. 3. There are three redundancies and the 
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; [X, 20 —134 —134]fm,]+ | 26,667 | , 

X;, 20 0 O}]] Mm, 0 
xX. 25 0 O}Lms 0 
Xa 15 10 0 17,500 
X, —5 20 0 10,000 

LX, 0 25 0 0 
x, 0 30 0 0 
Xp 0 20 10 10,000 
X, 0 0 30 OO 

where my M,, ms = M,, m, = M,, 


ind the division into groups is indicated by square brackets. This problem 
proved extremely difficult to solve with the iterative programme described 
in § 4, for reasons mentioned in that section. Using the second programme, 
however, no difficulty was encountered. 

In the case where the initial trial point was chosen to be the origin the 
results printed by the computer were as follows: 


rial First Vintmum 





f1on erlex ertex 
26,667 10,000 5,155 
10,000 5,155 
12,500 6,452 
co 22,5C0 15,278 
7,500 g25 
6,250 15,275 
7»500 19,333 
00 7.506 3,889 
7,500 15,333 
$0,000 35,790 
\ \ 
{ la Vp a b 
\ 
n=, Yq = Xy 
| \ . 
X, = X; X= X; 


It wili be noticed that movement from the first to the minimum vertex 
takes place along an ‘implied line’, as described in § 5.2.2(b). The time 
taken by the machine to produce the above solution was 40 seconds, of 
which approximately 12 seconds was devoted to actual calculation, the 
remainder being spent on input and printing. 
Example 3 

Che largest frame so far treated by the programme is shown in Fig. 4. 
rhere are 30 moment points and 12 redundancies. The programme success- 
fully found the same value for the minimum weight from a wide range of 


initial trial points, although it was found that there were several vertices 
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ossessing the same (minimum) weight. Starting from the origin m. = 0 
> ta) eS eS 1 


the solution time was 7} minutes, this being divided as follows: 


Input and initial printing : , : . 45sec. 
Calculation of first vertex : ; . 38min. 5 
Selection of a minimum vertex . oS » 25 
Printing of solution 20 , 
Se 








With the above starting-point 14 vertices were tested before the minimum 
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was attained. The average time of 15 seconds per vertex was divided as 
follows: 


Selection of new line 7, > a ' 3 sec. 
Caleulation of 8. : : ; ; : 2 
Matrix alteration . . , - : 10 


ia 
One interesting feature of this problem was that during the solution several 
cases of ‘coincident vertices’ appeared. At such a point more than r 
equalities are satisfied, but this does not affect the operation of the pro- 
gramme. The vertex is treated as two separate points joined by a line 
which, although of zero length, has a definite gradient associated with it, 
and this line is tested in the normal way. It is also interesting to note that 
an elastic analysis of the same frame took the computer approximately the 
same time (9). 


8. The calculation of collapse loads 

So far the problem considered has been the choice of a suitable set of 
fully plastic moments to support a given system of loading. The analysis 
which has been developed is, however, equally suitable for solving the 
related problem, in which the fully plastic moments are given and it is 
necessary to find the factor A by which the loading must be increased to 
sause collapse. From the linear nature of the equations of equilibrium it 
is apparent that the latter problem is equivalent to finding the largest factor 
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\ by which the known fully plastic moments of a structure can be divided 
before the structure collapses under a given load. 

To solve this problem by the analysis developed in § 5 it is necessary 
to re-define the variables X,. These will now be used to represent the 
quantities .V,/.%,, where ./, is the fully plastic moment of the member at 
the point « and is assumed to be known. Since the factor A will affect each 
point equally, it is natural to include all the points a in the same group. 
The symbol 4 thus indicates simply the value of « for which |X ,| is greatest. 
Dividing the equations represented by (1) by the appropriate values of @, 
gives a set of equations which may be represented by (2), provided that the 
coefficients A,; are suitably reinterpreted. 

With these slight modifications in notation the analysis proceeds exactly 
as before. Since there is now only one group, the function G becomes 


. My 

@ = (XZ; x 

M. 
That is to say, G represents the largest value of the ratio |M,\/.@,. For any 
given moment system, therefore, the largest factor A by which the fully 
plastic moments .@, can be divided without violating the condition 


M |r VW, is given by 


Clearly, the maximum possible value of A will be attained when the function 
(is minimized. As before, this minimum will normally occur at a vertex, 
where r equalities of the type |X, X;| are satisfied. At such a vertex 
there will be (r+ 1) values of X , of equal maximum modulus, corresponding 
to the (r+ 1) plastic hinges which are in general necessary for collapse. 
However, if A,,,,, is attained at more than one vertex, a collapse mode will 
exist which has less than (r+-1) plastic hinges. This will be associated with 
local failure of a certain part of the structure. 

\lthough the method described will determine the correct load factor 
ind final collapse mode in all cases, the path traced out in m,-space will not 
in general correspond to the true behaviour of the moments in the structure 
as the plastic hinges develop. An investigation of this behaviour requires 
knowledge of the elastic characteristics of the structure, as well as the 


equations of equilibrium 


9. Conclusions 
he analytical method described in this paper gives an exact solution of 
the general mathematical problem considered, while the digital computer 


programme enables rapid solutions to be obtained in individual cases. It 
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is of course true that there are many factors which affect the cost of a 
structure besides its weight, and in a practical design several different 
loading systems must often be considered. The computer programme, 
however, provides a method of solution so rapid and automatic that it may 
be of value in giving the engineer a rough guide in the initial stages of his 
design work. 
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SUMMARY 

A generalization of the incremental strain law given by Hodge, Prager, and 
Drucker is developed. The new law affords the possibility of several different types 
of loading regions in contrast to the single loading or unloading region of the older 
theories. The stress-strain relations can be explicitly inverted and minimum 
principles for the stress rates and strain rates proved in the manner of Greenberg. 
An absolute minimum principle corresponding to the differential equations of 
structural stability is developed. A particular evaluation of the coefficients indicates 
that the new law may have broader applications than the previous theory. 


1, Introduction 

Serious discrepancies between theory and experiment found in the appli- 
cation to buckling problems of the incremental strain or plastic flow law 
of Prager-Hodge (1) and Drucker (2) have been noted previously (see, for 
example, Pearson (3) and Bijlaard (4)). Some of the criticism has been 
directed towards the stress-strain law itself, although Onat and Drucker (5) 
maintain that small initial imperfections account for the differences. 
However, certain features in these previous results lead one to inquire 
whether more general flow or incremental strain laws are possible. 

The general form of the Hodge—Prager—Drucker law is 


; , of of. 4 
€ij H kt SKIT G . f , J Or» i f Ox > 0 
, 005; COjpy C0}, 
of . 
Ah 44 Fy: = Fn S 9, (1.1) 
sie 


where o;; is the stress tensor, ¢;; is the strain tensor, the dot represents 
differentiation with respect to a time-like parameter, f is the loading 
function ({é@f/éc,,\é,, > 0 representing loading), G is a scalar which may 


be a function of stress and strain, and H,,,, are the coefficients in Hooke’s 


+ Work on this problem has been partially supported by the Wright Air Development 
Center, U.S. Air Force. The main results were presented under a similar title at the 4th 
Symposium on Plasticity, Brown University, 1-3 September 1953. 
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law. Ifthe basic state of stress is pure compression (as in the plate buckling 
problem), the factors G éf/éc;; can be completely represented in terms of 
the quantity » = E,/H, the ratio of Young’s modulus to the tangent 
modulus in pure compression. This dependence upon a single modulus 
seems to be in part responsible for some of the difficulties. 

This paper shows that the above law can be generalized, within a similar 
framework of postulates, to one containing more parameters. It remains 
to be seen whether this new law yields any improvement in the buckling 
problem. 

Section 2 contains a development of the stress-strain law. It is first 
shown that the assumption of linearity of the strain rates as functions of 
the stress rates can be replaced by equivalent postulates which are more 
direct expressions of the physical nature of the material. The Drucker 
work-hardening criterion then permits an explicit representation of the 
coefficient functions appearing in the expression for the inelastic strain 
increments. The relation between this representation and the loading 
criterion is then examined. In the following sections, minimum and 
uniqueness theorems on the stress and strain rates are established and a 
minimum principle which characterizes the buckling load is derived. The 
non-elastic coefficients are evaluated explicitly for simple states ‘of stress 
in § 5. 

It should be noted that the stress-strain law to be established here has 
similarities to recent work of Koiter (6), although the derivation is based 
on a somewhat different point of view. The final forms are sufficiently 
different, however, to allow for a more direct proof of the minimum 
principles. 


2. Development of the stress-strain law 
The stress-strain relation to be studied describes the dependence of the 


strain rates €;; on the stresses, strains, and stress rates, namely, 


: ‘ - 
€44 €ij(Ox2> €mn> Fpq): (2.1) 


Let us assume that our work-hardening material is such that the strain 
rates have the following properties as functions of the stress rates: (a) the 
strain rates have continuous first derivatives as functions of the stress 
rates, (b) ¢;; = 0 implies that €;,; = 0, and (c) a change in the time scale 
does not affect the stress-strain relation. 

Assumption (c) implies that if 7 = At, then 


do,,,/dr). 


€i(Ox, € 


; mn? pa 
dr 
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kling Rewriting this relation in terms of the variable ¢ and using (2.1) we have 
ns of ] de;; : , a. , 99 
ee iv €ij| Fx» mn» om y Stal Ft» Emu» Fpq)- (2.2) 
lulus . i 
Equation (2.2) states that the strain rates are homogeneous functions of 
the stress rates of order one. Consequently, @é€;;/¢c,, are homogeneous 
. ‘ functions of the stress rates of order zero. However, a continuous, homo- 
rt geneous function of order zero must be constant. Therefore, the strain 
— rates are linear functions of the stress rates; and assumption (c) implies 
first that €ij Bi jx0 Fy (2.3) 
IS Ol where the B,,,. may be functions of stress and strain but not of the stress 
nore | rates. | . 
cker The assumption on continuity (a) can be expressed in a somewhat more 
the ? 
A 
rain OC 
ling | 
and | 
id a 
The 
ress 
| 
has 
~ | Gi 
itl Fig. 1. 
um 
useful form in the following manner. The new requirements can be most 
readily understood by using the geometrical language of stress space (see 
Fig. | for a two-dimensional sketch). The state of stress at a generic point 
ii in the material body is represented by a point in a six-dimensional space 
with coordinates o;,. Additional stress rates, starting from this state of 
stress, are represented by vectors with components ¢,, with initial points 
1) 3 at the stress point o;,. In Fig. 1 such a stress-rate vector is shown by the 
sin dashed line. We shall require that the €;; be continuous functions of 6,,, 
i but that their first derivatives with respect to ¢,, need only be continuous 
sis in regions which are convex with respect to rays emanating from the stress 
i. point (such regions are shown as I, II, III, 1V in Fig. 1). We further require 


that in any such region, the limiting values of the derivatives shall be 
uniquely defined as the stress-rate vector approaches zero provided that 
the path of approach is not tangential to the boundary of the region. Thus 
the strain rates, though continuous functions of the stress rates, have 
derivatives with respect to the stress rates which may exhibit jumps due 
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to the discontinuities in B,;,,in passing from one such region to a contiguous 
one. 


If the further assumption is made that Bjjjy = Hjjy+Djj, where H,;,, 
are the elastic Hooke’s law coefficients, then the D,;,; represent the non- 
elastic coefficients (D,;,,; = 0 in the unloading or elastic regions of the 
body). It remains to determine the form of the D,;,, for loading. 

From the symmetry of the stress and strain tensors, the coefficients Diy, 
have the symmetry properties D,.; = Dix, = Dz. We shall also assume 
the double symmetry property D,;,, = D,;;.. This assumption implies that 
only those coefficients are retained which contribute to the ‘increment in 
plastic work’, dW = dej;do,;; = D,;,,do,,do,;. Here the increment or 
differential do;; is defined as do;; = 6;,dt. 

If the Drucker work-hardening criterion (2, p. 412) is imposed, viz. 
€;;5;; > 9, with equality holding if, and only if, the plastic strain rates 


és 0, then — 
D449 F 55 Fj > 9. 


Thus D,;,,6;;6,, 18 a positive definite real symmetric quadratic form of 
rank six since but six of the 6;; are independent. The formal consideration 
of this as a form c,;7,2; > 0, ¢;; = c¢,;, leads to the conclusions that the six 
eigenvalues A”) of the form are real and non-negative. Furthermore, the 
matrix c;; may be written according to Wedderburn (7) as the sum of 
outer products of its eigenvectors b\”): 


6 
— (n) (nm) (m) f(n) (n) = ° ai J 
C= 2, by”, Mb) =X”) (m=n); =O (mn). 
at 


If we identify the x; with the ¢;;, ¢;; with D,;,;, and the eigenvectors b{” with 


‘eigenvectors’ A\”), A‘) — A”), in the same fashion, viz. if 24 = 64) = 6», 


then yy = 4Djy2 = Dyoye+ Dyyy2+ D221 + Day, and by” = AY) = Af? 
then D,;,, can be expressed as 
6 
Diya = 2 APAP, AP= AP, | ; 
n=1 (2.4) 


Aim Aim) =k (m=n); =0 (m #~ n)) 


We shall speak of the A{? as the eigenvectors of the matrix D,,x). 

Equation (2.4) cam also be given a simple geometrical interpretation if 
A} is considered to be the normal to a surface element in stress space. 
At a given point in this space, a corner is formed by the intersection of 
these elements, and a set of regions is formed which are bounded by 
consecutive elements. The number of such regions depends on the number 
of non-vanishing eigenva' 


Various states of loading are possible depending on the number of terms 
> Lo) 
(n) (n) « “7 j k adi iti nll anecify whic 
A”, Aj}? appearing in Dy 54g: The loading conditions will specify which of 


— 
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the terms are to be retained. Thus six loading conditions will in general 
be required. We shall assume that these are linear in the stress rates and 


are of the following form: 


By? 6, 0 implies A}}? is present, 
By? 6, 0 implies A} is not present. 


Consequently, if By c,, < 0 for all n, the stress-strain law is purely elastic. 
Such will be the case for complete unloading. We shall also postulate a 
continuity property which is a generalization of that previously used by 
Handelman, Lin, and Prager (8) to describe behaviour under neutral 
changes of stress. We shall assume that if the stress rates ¢,, are such that 
3”) G approaches zero through positive values while the signs of all other 
\") ¢. remain the same, then A‘? Aj 6,, approaches zero. Consequently, 
if BY 6. = 0, then A}? a, = 0. Thus Aj? and By” are scalar multiples 
of one another. Since A}¥” is determined ‘aly to within an algebraic sign, 
we may replace the loading condition BY” 6,, > 0 by Ajj 6,4 > 0. The 
final stress-strain law can be written as 


: . (n) 4 (n) 
€ij Aj 543 Fy 1. > A‘ ; A ; ? ou | 


nm 
a oe ee 
itki E, ik°jl E, ij “kl 
where only those terms are retained for which Ajj 6,, > 0. Here E, is 
Young’s modulus and v is Poisson’s ratio. Equation (2.5) may also be 


written as 6 


. ° ( ( . (n) 
€; ii 53 4g +4 & A ij | AX? Gy T Aj j.\)- 


v tj Fe 

The stress-strain law (1.1), considered in the previous section, is seen to 
be a special case of (2.5) with A} VG of 00;; and all other Ay 0. 

Essentially, the stress-strain relation (2.5) provides a ‘corner’ at every 
point in stress space, where not one loading surface is assumed to exist now, 
but where six surfaces intersect orthogonally (see (2.4)). Indeed, at each 
point six surface elements, not necessarily integrable, need only exist. The 
stress-rate vector at each point may be oriented into a total unloading 
region or into a region where all, one, or any combination of the six quan- 
tities A\") é,, may be positive. A two-dimensional sketch of such a corner 
is shown in Fig. 2. For the sake of simplicity only two terms, Aj}, A\?, 
are taken in the stress-strain law. The surface elements and their corre- 
sponding normals are indicated in the drawing together with the various 
possible loading and unloading conditions. We have also introduced the 
notation L, to refer to a region for which A}} é,, > 0 and U, to refer to a 


region for which A}} 6;, , with corresponding meanings for L, and U, 
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Our previous discussion shows that if the surface elements are actually 
integrable, then the inelastic strain rates are found as a linear combination 
of gradients of potentials. Furthermore, these potential surfaces corre. 
spond to the loading surfaces. The physical significance of the integrability 
requirements is being examined together with possible examples of the 


types of loading conditions to be expected. 
Unlike Koiter’s (6) stress-strain law, (2.5) can be inverted, if the assump- 
tion of plastic incompressibility is imposed in each possible loading region, 











4 
C2 
On 
Fic. 2. Region L, L,: Af}6;; 0, Af?) ;; 0: 
region L, U,: Afj)6;; > 0, A?6;, < 0; 
region U, U,: Ajj)6;; < 0, A? 6;; << 0; 
region U, L,: Ajj)é;; << 0, A6;; > 0. 
rh (n) § 
viz. Ay 3,; = 0, n | 


oa 6. Because of the ort hogonality,+ Greenberg’s 
(9) method may be applied to yield the inverted forms 


6 
: 5 2 .: (n) R(n) = 
Oy; Al yin € ji é, A ij By; Ek 


n 
E v 
1 40 Ps ~ o ~ 
Ain | in Or; O51 + Dp oti On, ’ 
‘ — = 
BM” Ey {l+1 4(n) 4(n) 4(n) (2.6) 
kl l i = | E T 4 pa 4 pa 4 kl . me 
40) 


with A\Pe,;, > 0 or BY é,, > 0 as the loading criterion. We note that 
(n) Rin) : (n) A (n) 
A\! By; By 44 je] - 


3. Minimum principles for the stress and strain rates 
Minimum principles of the generalized Hodge—Prager type for the stress 


and strain-rate boundary value problems hold for the new stress-strain 
+ The authors would like to express their appreciation to Dr. Bernard Budiansky for 


pointing out the importance of orthogonality in allowing for an explicit inversion of the 
stress-strain laws. 


He has noted that the difficulty of inversion is present not only in 


Koiter’s relations but also in the slip theory relations proposed by himself and Dr. Samuel 
Batdorf. 
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relation. They are of the type given by Greenberg (9). In our discussion 
of these principles, we shall only give those portions of the proof which 
differ essentially from his. We shall restrict ourselves, however, to a single 
loading. Koiter (6) has also proved similar principles by methods modelled 
after Hill (10). The chief difference in what follows lies in the fact that 
Koiter’s law cannot be inverted explicitly and hence his strain-rate principle 
is not amenable to Greenberg’s method of proof. A by-product of the 
computation is an inequality concerning the coefficients of the stress-strain 
law and admissible strain rates. 

The problem which we shall consider here is the following. Let us assume 
that the stresses o,; are known throughout the volume V of a body. On a 
part S, of the surface, the displacement rates w%; are prescribed; on the 


remainder S, of the surface, the traction rates T; = 6;;n; are prescribed. 


It is required to determine the 6,; or the €;; throughout V. 

An admissible set of stress rates o*; satisfies the equilibrium conditions 
in V and the stress boundary conditions on S,. The corresponding strain 
rates €*, are defined through the stress-strain law (2.5). The function J(6%,) 
is defined as 

ii €ij T¥u, ds. 


j Si 


I(6%) = 4 | oh & dV 





An absolute minimum principle for the stress rates can then be stated as 


foll WS 


STRESS-RATE PrincrpLE. Among all admissible stress rates o},, the true 
solution G;; to the boundary value problem is that set which minimizes I(67;), 
and the &;; are unique ; that is, I(6%,)—I(6,;) > 0, with equality holding if and 


only if G:; Oj; 


The proof of this principle will not be given here since it follows directly 
along the lines of the method of proof given by Greenberg. It is only 
necessary to observe that the integrals in his proof are replaced by finite 
sums of integrals, one for each type of loading region. 

For the strain-rate problem, an admissible set of strain rates ¢*, are any 
derived from velocities u* satisfying the boundary condition on S,. The 
corresponding stress rates 67, are obtained through the inverted stress 
strain relation (2.6). The function J(é%) is defined as 

J(c%) = 3 | che dV— | T,at dS. 
i 82 


We then have the following 


STRAIN-RATE PRINCIPLE. Among all admissible strain rates €},, that set 


€;; Which solves the boundary value problem is the one which makes J(€3;) an 
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absolute minimum; that is, J(é};)—J(€;;) > 0 with equality holding if and 
only if &; = €;;. 

Let €;; = efj—é€,; and u, = u*—1,. By a process exactly similar to that 
used by Greenberg, with the addition of the interchange of sum and integra] 
over each of the loading regions, it can be shown that 


- 6 

ok . . 5 * rs 

J(%)—J(E,;) > $ | > AMBP & dV + 
tbe*=* 

6 


a 22 4 + 2 , 
T } Aina €i3 © — 2 AW”) BYD E,5 Ey d\ . 
} | ; | 


A n 
J 


The notation L, U refers to loading or unloading regions corresponding 
to €,;; L*, U* refer to similar regions for <7. A double suffix refers to a 
region where both relations hold. The first integral on the right is always 
positive since A{¢* and Bie*, have the same sign. The use of generalized 
inequalities of the Cauchy type is needed to show the positiveness of the 
second integral. Assuming for the moment that the integral has been 
shown to be positive definite, the argument is easy to complete; for, by 
this assumption J (€¥;)—J(€,;) > 0, and, if ¢% = €,,, J(é%) = J(é,;). If the 
latter holds, then the second integral on the right must vanish. By its 
positive definiteness, it will vanish only if €,; = €4;—é€,; = 0 throughout J’. 
This completes the argument except for the proof of positiveness. Since 


H 3 = = E, 2 2 ' Vv = + 
ijkl ©ij Ek lt p\ ST 7. Dp <4 45 


and Ey vé;; €;; (1-+v)(1—2v) > 0, E,/(1+v) > 0 


~€ 


the problem reduces to showing that 


+ 6 
—{/jJtyp 1 
= = | (n) 4(n) (n) A(m)= = ~ 
| (e,2— S BE --— Any AS A ij Aji €5j Ex dl 0. 
ny “0 

The integrand will be shown to be positive definite. For notational pur- 
poses, let (1+-v)/EF, x > 0 and let the €,; and A‘ be replaced by the 
nine component vectors ¢,, and A\” in the same way, viz. if €,, = €,, then 


AS) = A™. The inequality to be proved is thus 
9 a 6 9 (n)? 1 9 \2) 
> &— > {lot ¥ AP] "(> AP) "| > 0. (3.1) 
k=1 n=1 k=1 k=1 


If (3.1) is cleared of fractions, we get 
6 9 2 6 6 9 1 / 9 ‘ 
of (n)?* 2 4 ia (p)? (n) oi 29 
I] jot > Ag) 2 &— 2 {TT (e+ & AP”) ( FAP a) SO, (3.2) 
n=1 k=1 k=1 n=1'‘\p=1 k=1 k=1 
a ; 
where | |’ indicates the product of all factors except for p = n. Expansion 
p=1 
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of this expression gives a polynomial of degree six in «; the proof of positive- 
ness consists in showing that the coefficient of each power of « is positive. 

Suppose that the vectors aj) of m components each, i = 1...., n, 
k = 1...., m, are not linearly dependent. Hardy, Littlewood, and Polya (11) 
state as the generalized Cauchy inequality of order n the following: 


m m m 
(1 l * af 8) lS (1) ,(n) 
> aPaP Yaa . . . Yapap | 
k=1 k=1 k=1 
m man 
‘ (2) ,(1 4 (2) ,,(2) 
ae aj” at: > a)” a}; 
| | >0. (3.3) 
m m 
N ama 1 ; : : 5  y aw a\™) 
me \ A om ‘. 
KI KI 
In the problem considered here, n 7 and m = 9, the vectors a? being 
BP 0s A$), and e,. Moreover, since Af”) A” = 0, m 4 n, it is also true 
that . 
> AMAM — 0, msn. 
k=1 
If a}! e, is always kept in the above inequality, but any number of the 
Ay”), m aa 6 are used, various inequalities of the type (3.3) of order 


up to seven can be written. All of these of concern here have zero elements 
except in the last column, last row, and main diagonal. If the Ay” are 
considered as ordered, then the number of different inequalities of order 


(q+1), q = 9...., 6, is 6!/[q!(6—q)!]. 


Consider a Cauchy inequality of order (¢+1), 1 <q < 6: 
% (n)? An) 
> Aj 0 . - : > Ay” «, 
k: k 
Zz Aj" 3 
. a a ' > 0, 
0 . . DA’ > Ape, 
k ke 
. f(m pa 4 (nq) , 
ye : . saR™ &j €k 
; f Lz 
with ny < m, <... <n, being q integers selected from the set (1, 2,..., 6). 
Evaluation of the determinant gives 
Nq itg ig 9 
‘ (n)* Se % , Ap)? An) i . 9 
IT > A; > «i > (II > AF (> Av €) 9, (3.4) 
r i kK k N=ny prHny, iC kc 


where the prime again denotes the absence of the term p = n from the 
product. An inequality of the form (3.4) exists for each of the 6!/[q!(6—q)!] 
orderings of the integers (1,..., 6) taken q at a time. 
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Denote by C(n,7r)D™ the sum of all possible ordered products taken r 
at a time of the D™, m = 1.,..., n, and by C,(n—1,r)D™ the sum of all 
possible ordered products taken r at a time of the D™, m = 1.,..., (p—1), 
(p+1),..., n. If all inequalities of the type (3.4) are added for each q, then 


0(6,q) > Ag") > <j. In the terms 
£ 


k 





the term involving > «7 becomes 
3 


9 
involving (> Aj” «,) , each such term for m = 1.,..., 6 will appear once for 


each ordering in which the integer m is included. It will be multiplied by 


the product of the other (q—1) ¥ Av’, n 4 m, appearing in that ordering, 
E ; 


and hence in the total sum it will be multiplied by the sum of all possible 
ordered products of the other five p32 Aj”, n ~ m, taken (q—1) at a time. 


Hence the summed Cauchy inequalities of order (¢+1), 1 <q < 6, can 
be written as 


(m)2 2 . Y (nm 2 (n) \ » a 
[C(6,9) ¥ Ayn") ¥ > [6.15.4 1S Ag |(E Ay e,) > 0. (3.5) 


If the expression in (3.2) is expanded, it is found that the polynomial 
can be written as 


x8 Y ef. >) "| 


k r=0 





C(6,6—r) > Ay| > i-— 
k 


6 9 
— > [e(6,5—r) ¥ Ag’|(¥ Ale.) I. 
n=1 k k 


Since ¥ «7 > 0, and since the coefficient of «” is seen to be the summed 
k 





Cauchy inequalities (3.5) of order q+-1 — 7—r, it is seen that each coeff- 
cient is positive and inequality (3.2) is true, thus completing the proof of 
the strain-rate principle. 


4. Minimum principle for structural stability 

It has been noted by Prager (12) that the differential equations for the 
general problem of structural stability also apply beyond the elastic limit 
for a wide class of stress-strain laws. The law proposed here falls within 
these requirements. In addition, Prager developed a corresponding varia- 
tional principle which has been extended by Weiss and Handelman (13) in 
the elastic case to an absolute minimum principle. It will be shown here 
that this absolute minimum principle applies equally well to the stress- 
strain law (2.5) and hence to the Prager—Drucker law (1.1). 

It should be noted that the stability problem is considered here only in 
the classical sense; that is, the Shanley hypothesis is not used nor is the 
effect of initial imperfections taken into account. Specifically, we consider 
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an initial equilibrium state of stress o,; and ask for the value of the multi- 


plier A such that Ao,; is a stress system in indifferent equilibrium. By this 


we mean that for the state Ao,,;, there is a non-vanishing set of perturbation 
stress rates ¢;; and corresponding strain and displacement rates, €,; and w,, 
which also form an equilibrium system. 

Only that portion of the proof will be given here which is essentially 
different from that of the elastic case. The remainder of the argument will 
then follow as in the previous result. The minimum principle states that 
the multiplier A is given by 


A min {D(u)/H(u)', 


u 


where . Z 
D(u) | Cysnr ss Ex OI | Crna Mig Meg EV, 
" Vv 
H(u) | o;;(€;.; € Up g Up) AV 
j 
| | o (3%, Uy j —_ Usk Ui Wy 4 Usk = Us tty, ;) dV. 
; 
and wi, is an admissible displacement rate. Let 
Diu. v) | O54 5 Pag di 
v 
and : 
H(u,v) 4 | o LOUK Ung — Vik Uj,k Uni Pinging leg) a A 
Then 
D(u-+-v) D(u)+-2D(u,v)+D(v) and H(u+v) H(u)+2H(u,v)+HA(v). 
Let wu i, be the true minimizing function so that D(w) = AH(u) and 
w= w u-t+ev = u,+e#,; be any trial function, with e an arbitrary scalar 


multiplier and #, an arbitrary, permissible displacement rate vector. Since 
the loading regions for uw and w will not in general correspond, let integrals 
involving w loading regions be denoted by an asterisk. Then 

D*(w) > AH*(w), 
or D*(u)+2eD*(u, v)+eD*(v) > ALH*(u)+ 2eH*(u, v)+eH*(v)]. (4.1) 


Since the H functions do not involve the stress-strain law, H* H. Now 


kl i kl 


* - 6 
D*(u) = | Hightigs ting V — |X AP BY ti, stig, AV 
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where the integrands have been suppressed for brevity and L,, 
to loading or ines regions for the nth term A’. It is easy to show 


0 < AW BYP ti, 5) < —EAP BY ai, ; 


rT <—e co A\™ BY D ti; 5 


D*(u) < D(u)—e | b> AM BO a; ,6,, dV. 
CU 1 


Substituting in (4.1) and using the relation D(u) = AH(u), we get 


> Al 2eH(u, v)+H(v)]. 
Similar operations in turn on D*(u,v) and D*(v) give 


- 6 
2«D(u,v)+eD(v)+e | Y AMBY 
7 1 


To compare these quantities, we must first eliminate the dependence on 
e of the regions U L* and LU* and, secondly, must express D(v) in terms 
of the v loading regions to compare it with H(v). Since the integrals over 
UL* and LU* are integrals of sums of squares (see (2.6)), the regions may 
be increased to the whole volume and the inequality can be only improved. 
If two asterisks denote the loading region for v, 


- 6 » 6 
D(v) = D**(v)— > AP BYP &; ; b,,dV | > A™ BY 6, 56 
“ n=] . 1 


Thus the left-hand side of inequality (4.2) is less than 


“| >> A) BY 6; ; 


| > Al”) Bo, 


Since D**(v) > 0, and since all integrands in 


the last three integrals are positive, the inequality (4.2) can be written as 


—AH (u,v) and B > 


The quantities A and B do not 
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depend on ¢ but only on w and the arbitrary vector v. The only way that 
(4.3) can hold for all « is for A to vanish; that is, 


D(u,v)—AH (u,v) = 0. (4.4) 
Once (4.4) has been established, the remainder of the proof of the minimum 
principle is exactly the same as that given by Weiss and Handelman and 


will not be repeated here. 


5. A particular evaluation of the A‘? 

As an example, we show that the proposed stress-strain law has a wider 
range of possibilities than that given by Prager and Drucker by evaluating 
the non-elastic coefficients for a state of simple compression, viz. 

O11 —G, < 9, 
all other o,; = 0. 

Consider various possible stress-rates starting from this initial state. In 
the case of additional simple compression, 6,, < 0 and all other 6,; = 0. 
If E denotes the tangent modulus in compression, E, Young’s modulus, 


and », E/E, we find 
Go E 6 -1 
11 i 0 7 7 (n) n) 
i B Eyj1+Ey y AQ AY)". 
€11 1 n=1 
6 
T 4 (n) (n y 
Thus > Ay’ Ai? (n,—1)/Ep. 
n 1 


If one now assumes that there are no permanent volume changes, and that 
compression in the x,-direction leaves the material isotropic with respect 
to behaviour in the 2,- and 2,-directions and produces no shear strain 


increments, then 


8 6 6 
n) (n) é . (n (n » ( (n) 7 
pI Ay? Ay ~ 2, Ay’ Aji’ = 2, Ay} Ay; (m,—1)/K, 
1 n n 
. 6 6 
. n) A(n) , «6A(n) Atm Ain) A(n) , at 
> AW AL 2 Ai AS) 2 AG AW=0 (t #}). 
n=1 n n 


Let us now consider the stress-rate tensor, 6555 for a general plane state; 
that is, ¢;, = 0, but otherwise arbitrary. Utilizing the results previously 


obtained for A\™, we find 


6 
é19/€2 = Eyl 1+v+2E, > Agy A\p] 
n=1 


-1 


(hus a shear modulus, G,, = 6 9/€,9, arising in a state of initial compression 
can be introduced. Further, if we let n,. = E,/G,., then 


, ] 
> Aj? Ay? op (ne 1—v). 
n=1 assy 
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Unlike the Prager—Drucker law, this generalization allows for a tangent 
shear modulus other than the elastic one. Furthermore, the loading con- 
ditions are sufficiently flexible to allow for an elastic shear modulus if only 
shear stress rates are present, where: s a plastic shear modulus can occur 
if a suitable combination of shear stress rates and compressive stress rates 
exists. Such a situation would require only that 


AY) 642 <0 for all n, 


whereas AM 64, + AD Goo+2A%) Gy. > 0 
for some n. 

Proceeding in the manner outlined above, one can evaluate the basic 
sums of the A} in terms of other parameters. It appears that there are 
six such parameters. 


6. Concluding remarks 

The results presented here are an attempt to show that the general 
physical postulates for a work-hardening material imply a more general 
stress-strain law than that given by Prager and Drucker. In some respects. 
then, they should be regarded as highly speculative. Nevertheless, there 
are several features which might be worth noting. The relation between 
linearity and independence of the time-scale has been clarified somewhat. 
The proofs of the minimum principles show that the mechanism responsible 
for their validity is broader than the Prager—Drucker stress-strain law. 
Finally, the minimum principle which describes stability in the classical 
sense has been extended beyond the elastic case. 
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THE DEFLEXION OF A SPRUNG RAIL 
By H. L. COX (National Physical Laboratory, Teddington, Middlesex) 


[Received 12 September 1955] 


SUMMARY 
Formulae are derived for the deflexion under a concentrated load of a flexible rail 
supported on sprung chairs uniformly spaced. 


1. Introduction 

THE deflexion under concentrated load of a flexible rail supported on sprung 
chairs has previously been treated by Timoshenko (1) and Inglis (2). 
Inglis’s analysis is based on the assumption that the rail is supported on 
eight chairs only, four on either side of the concentrated load, the small 
fraction of the wheel-load which is transmitted to chairs yet more remote 
being disregarded as negligible in magnitude. Fora very stiff rail on rather 
flexible chairs this approximation may prove insufficiently accurate. On 
the other hand, Timoshenko’s analysis relates to a rail on a continuous 
elastic foundation; and this representation, which is certainly adequate for 
a stiff rail on very flexible chairs, may be insufficiently precise for a flexible 
rail on rather stiff chairs. 

The present note describes a complete solution for a rail supported on 
any number of chairs, from which the solution for an infinitely long rail is 
inferred. The effects of shear deflexion of the rail and of rotational restraint 
at the chairs are disregarded; but these influences could readily be repre- 
sented by slight modification of the basic analysis. 


2. Statement of the problem and notation 

The rail having moment of inertia of section J is supported on N+] 
chairs spaced h apart, so that the total length of the rail 1] = Nh. The 
defiexion under load W midway between two chairs at which it might be 
simply supported would then be 5, = Wh?/48EJ. The deflexion of a single 
chair carrying by itself the same load W is denoted by 6,. In all that follows 
the characteristics of rail and chair are represented completely by the 
deflexions 6, and 6, under unit load. 

Distance along the rail from one end is denoted by x and the deflexion 
at x by y. The latter is represented by the form 

yY = ¥ +o(a—$l) + > A, sin(nzx/l).F (1) 
; 

+ Formally the constant term y, and the si term % may themselves be represented 

in sine series; but the explicit form is preferable, because it represents the differentials of y 


in more convenient forms. In fact A, as written converge as 1/n*, so that the shear load 
proportional to d*y/dz* is still represented by a convergent series. 


(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 3 (1956)] 
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3. Solution of the problem 


The total energy of deformation is 
l ' 
. . od ats So N , cs 
SEI | (d*y/dx*)? dx +3(W/8,) > y2—Wy,, (2) 
0 r=0 
where y, is the value of y at the rth chair (x = rh) and y, is the value of y 
under the load W applied at x = a. By substitution of the form (1) assumed 


for y, the total energy for unit load (W = 1) becomes 
of - » 2 
—— » ntd42?+— S [y.+yl{(r/N)—}}- +> A, sin(nzr, y)}* 
192A OR Land «0, ind 
~ on r=0 


—{y,+y(a—}l) + > A, sin(nza/l)}. (3) 


By differentiation with respect to y, and equating to zero we find 


(N+1)y.+ > A, cot(nm/2N) = 6,. (4) 


Cc 
n odd 
By differentiation with respect to % and equating to zero we find 


{(N+1)(N+2)/12N}bl—4 >} A, cot(nz/2N) = 8,(a—4l)/l. (5) 


neven 


By differentiation with respect to A, and equating to zero we find 


4 . r=N 
TT Op . . yr 
a) ee i, Y. sin( +-pl{(r/N)—4}sin(nar/N) + 
O6N 3 0, kan _ 
r=0 


> A,, sin(mar/N)sin(nar N)| = 6p sin(n7a/l) 


m 


or, performing the r-summations, 
(74/96.N*)n4A,, + (5;/8.){y, cot(nm/2N)+4N56,} = dpsin(nza/l) (6) 


for n odd, and similarly for n even with —43yl for y,, where 


8, = A,—Aoy-ntAenin—Aan-ntAan-n—-o+s ete. (7) 
since 

LIN ifn = 2Np+m, where p is any integer in- 
r=) aa _ | cluding zero, and values of m and n mul- 
> sin i sin — 7 tiples of N are excluded ; 
tes ; | 0 ifn 4 2Np-+m, or if m and n are multiples 

of N (3). 

Since cot{(2rN +-n)(7/2N)} cot{7r+(na/2N)} = +cot(nz/2N) and 


Oonp+n On = —ava—n 
formula (6) may be summed to yield 


(74/96.N3)5,, + (8,/8,){y, cot(n7/2N)+4N5,}S, = 5p 8), (8) 


n> 
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for n odd, and similarly with — 31 for y, when n is even, where 


1 l ] m7! 2-4 cos(an N) 

S,, = a ee - — ——— -o. SS Oo 
nt  (2N—n)*  (2N-+n)! 48N4 sin4(zn/2N) ” 

yy sinnd sin(2N—n)d . sin(2N-+n \p “i 

A == « -_—+-— eees Cc. 

‘ n' (2N—n)! PNG ny 
and ¢ = (za/l). 
The summations represented in 6, = A, —Aogy_,+Aoy.,—-+) ete., are 


over all values of n except multiples of V. For the latter there is, of course, 
no contribution to the energy from the supports (y, = 0 for all n multiples 
of N), so that 


(74/96.N%)(pN)*A,y = 6,sin(pN¢) simply, (9) 
where p takes all positive integral values. 


Recalling that cot{(2rN-+n)(7/2N)} +cot(nz/2N), formulae (4) and 
(5) may be recast in the form 


(N+1)y.+ > 5, cot(nz/2N) = 5, (10) 
nodd 
and {(N+1)(N+2)/12N}bl—3 SY 58, cot(nz7/2N) = 6,(a—3fl)/l, (11) 
neven 
where the summations are now restricted to values of n < N. 


For given values of (a//), N, 5p, and 6,, values of 5,, follow from formula 
(8) and values of y,, ¥l, and A, (n 4 pN) follow from formulae (10), (11), 
and (6). The deflexions at the chairs are then completely defined, but to 
complete the deflexion of the rail between chairs the terms in A_,,, defined 
by formula (9) have also to be included. By expressing the contribution 


from these terms to the total deflexion in the form yy = > A,,y sin(pNzz2/l), 
p 


it will be found that y,, would result from the loading 
= EI (d*yy/dx*) = (2W/Nh) > sin(pN¢)sin(pzax/h). 
> 


Then by substituting a = (r+v)h, where 0 < v < 1, so that 
Nd Nzaa/l = zalh n(r+yv), 
it follows that wy = (2W/Nh) > (—1)?’sin(pzv)sin(pza/h). This is the 


D 
Fourier representation of loads W/.N acting up and down in adjacent bays 
at (r+v)h and (r+2—v)h. In effect y, is periodic but antisymmetric in 
adjacent bays; and in the loaded bay it gives the deflexion of the rail when 
simply supported at the two adjacent chairs and loaded by (1/.)th part 
of the unit load applied at the section vh. When N is large this contribution 
yx to the total deflexion becomes negligible. 
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4, Application to a rail on three chairs 
When V 2, the only value of n in formulae (8) is n 1; then S, = 74/96 


and 


sneod | a4 ~~ . 
Ss) $ (-1 = ¢ 3- 4? for —47 < ¢ < }z. 
odd p* 96 7 i 


Hence, it is convenient to substitute 2a/l 24/7) and maintain 
P 


1 when formula (8) reduces to 


p 


15, +(3 7/8. Ye+8,) = 43p p(3—p?). (12) 
Formula (10) yields 3y,+6, = 4. (13) 
and formula (11) dl 8,.(p—1). (14) 
Substituting J (6./5,,) for brevity, 
(5,/8,) 4(2—9p+ 3p?)/(3J + 16), (15) 
(y,/0,) (J +-8—12p+ 4p?) /(3J + 16). (16) 


Then from formula (6). 


T7680 . nmap 2J+24 8p . 
4 Rain r. “ ‘P * pP for n odd 
win! | 2 3J +16 
and from formula (9), ‘ (17) 
7680, . Nn7p ; 
A. Rsin— for n even 
mn ys 


The coefficients A, are, however, only needed to describe the deflexion 
of the rail between chairs, and it may sometimes suffice to restrict attention 


to the deflexion at the chairs only. The deflexion y, at the rth chair is 


; ] *. . ar 
Yy,. Y, wl )4 S 6,, sin 7? (18) 


n 


and hence for the case of three chairs 


4/5 y.— sul J +8—12p+4p3 : l—p (19) 
Jol 5 “3d +16 :° 
) ; ‘ | 9 § 3 
Y3/9, os 1 : x Pa (20) 
5 3J +16 
? Ly J+8—12p+4p? 1- 
Yo/O¢ Jet? be - ne c. (21) 


~ 


3J +16 
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The partition of the loads between the three chairs is indicated in Table 1, 


TABLE | 


Partition of load on three chairs 

















p |e t+ | 3 i 
Yold_ | 40f+256 | 34J+177 | 28J+104 | 22/+43 16] 
44/8, 16] | 16]/+94 16]+176 16]+234 | 16]+256 
42/5, —8] —2]—15 4J]—24 1of—21 | 16] 





All values are divided by 16(3J +16): when p > 1, substitute 2—p for p 
and interchange y, and yp». 


5. Rail supported on many chairs: general formula for deflexion 

The results derived in the previous section illustrate the complexity of 
the complete solution for any finite number of chairs. In such cases, if the 
deflexion at every section of the rail is required, the method of analysis 
used by Inglis (1) may indeed be preferable. However, if attention be 
restricted to the deflexions at the chairs, application of formulae (8), (10), 
and (11) suffice and the solution of simultaneous equations is avoided, 
Moreover, application of formulae (6) and (9) to determine the deflexion 
at specific points, as, for instance, under the loador loads, is straightforward. 
At the same time, the main purpose of the present analysis is to arrive at 
the solution for a rail supported on very many chairs, which solution is not 
accessible by the approach used by Inglis. 

When the rail is long, y, and % both tend to zero and formulae (10) and 
(11) are not needed. Then we may write the deflexion y, at b (= ph) due 
to unit load at a (= Ah) in the form 


Yo = > A,, sin(nur/N) (22) 
n=1 


and, since N is very large, we may write nz/N = @ and treat @ as a con- 
tinuous variable, so that 











yp, = (N/m) { A, sin pO dd, (23) 
0 

where, from formula (6), 
NA,, = 968 ,{sin AO—4N(5,,/5.)}/04, (24) 
and from formula (8), N58, = 96S} 5 ta-FP (25) 
where F = (1—cos@)?/4(2+cos 6) (26) 
; " ‘ ” sind@ sinA(27—8@) . sinA(27+-@) 99 
ee, « "a ae OO 


oe see 
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so that NA,, = (968 p/6*)[ sin A9—48S8)/{J +(1/F)}] (28) 
965; , - 48S) \sinué ; 
and Y,, oy | sin — sa 54 dé, (29) 
0 


where J = 6,/dp. 

Any integrand which includes a trigonometrical factor with argument 
(A-+-)@ oscillates with period 27/(A+.) in 6. When both A and yp tend to 
infinity these integrals are all zero. Therefore formula (29) reduces to 


(30) 


Ré : a > ZB f( 9 
m —e | enst—p0— 2. S ee. 
’ l+JF 7 (0-+-2n7)* 64 
0 aun 


where the summation is over all integral values of n both positive and 
negative and including zero. 
Now sa 
” N cos(na+y) ] & | ~ cos(na+y)) 
BY ~~ 6B 2. nth 
alln " ?) B \ alln ; B J 


and 


| Bin ' pon | 


cosy | S 26 COs Nx COS y_ c feenontonead 5 
B?—n? Gul pe—n? | 
n 


4 N=H y . 

cos y 2 7 C ~~ £COS8Na 
—— On y — Gn? - een @ 
8 Ox pa B?—n? 


n= 


) ) 


— it 8) 2) 
ally i 


XS cos(na+y) cosy | ~ (cos(na-+-y) , cos(na—y)) 
i 
y a) 
U 


1 


Moreoyv er, 


COS Tia 7 Cos B(z - a ] 
= : 7s... = for0<a< 27 
9 R nm” B sin 7B Bp? 
T 


and therefore 


* COS(?a--Yy TT 
7 y) — — cos{B(z 
<~ n+p sin 7B 
a 7m 


for all values of y. 
And finally 
N cos(na y) O° | 


(7 e) 


cos{B(7—a)+-y} | 





6 sin 7B 
7 (2(2+cos27B) 3(a/7)? (:) 
6 


= per —vit 


—_—_ 
alln 


|" sin’xB ——srsin°8 


? ‘ = ie (2) cot mB} sin 8—y)}. (31) 
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In applying this summation tc formula (30) we have B = 6@/2m and 
y = (A—p)@; but for ~« we must substitute 27v, where v is the excess of ) 
over the next lower integer. Then formula (30) may be written 





485, 
Y, = E | cos(A—.)0 — 
7 a 
0 
{(1—G,)cos(v—A-+-)6-+- H, sin 6 sin(v—A-+-.)@) | dé (32) 
| cen: << ae flor 
where G,, = {3—v(1—cos @)\v?(1—cos @)/(2+ cos @) 
and HH, = {3—v?(1—cos @)\v/(2+ cos 8). 


Each of the functions F,, G,, and H, is unaltered in value by substituting 
2n7+86 for 6. Therefore, if the whole range of integration from 0 to infinity 
be divided into subranges 0 to 7, 7 to 27, ete., formula (32) may be con- 
verted to the form 





4895p [¢- , 
Yo EY [fA—n) {UL —G,) fv A+) + SA fv—A+p— 1) 
0 
SH, f(v—A+p+1)}/(1+JF)] dé, (33) 
a , l * cos{2ane + (A—p)6} 
where f(A—p) 3 ya ee j 


and « is the positive fractional part of A—p. 

The forms for the other functions f( ) are similar, and thé fractional part 
v—e is the same for all three. Since A and p are always interchangeable, we 
may take A—p = M-+-e, where the integer M may be positive or negative, 
and assume v—e > 0. Then 

48 f(A—p) = {((1—G@,)cos M6— H, sin 6 sin M0} / F, 
48 f(v—A+-p) = {(1—G,_,)cos M0+-H,__ sin @sin M6} F, 
and 
48{ f(v—A+p—1)—f(v—A+p+1)} 
2 


F {—(1—G,_.)sin M6 sin 0+-H,_. cos M@ sin?6}. 
Then the formula for the deflexion at vh in one bay due to unit load at 
(v—e)h in the Mth bay in the direction of v decreasing is 
y ok | [{a—4@,)1+J F)—(1—4@, )(1—G,_,)— H,, H,_, sin?@}cos M6 

0 


dé 
G ‘sin 6sin Mé@ —s 
(1 »—e)H,\sin @ sin MOl iG TF) 


(34) 


—{H(1+J F)+(1—G,)H. 


v—€ 
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and 6. Particular cases 
; 

of ) (a) Loads on chairs 
; Since G, and H, are both zero, when v and « are both zero formula 


(34) reduces to . 
Yur ok {J /(1+J F)}cos M@ dé. (35) 


« 





’ 
29 Thus the proportion of unit load applied over a chair which is carried by 
the Mth chair on either side is 
Yas lf cos Mé dé a 
Pu — a ° (36) 
5 a J 1+J(1—cos 6)?/4(2+-cos @) 
Ing 0 
. ’ ° . . 
ut} \ll the integrals (36) may be evaluated in algebraic forms (e.g. by 
On- § substituting w tan (4@)), but the resulting forms such as 
fs) } ot. ») romp Q7 
Do = (3!+(1+-J)-}/[2(24 341+ J) HF (37) 
are complex, and for numerical evaluation direct integration of formula 
35) is preferable. Values of p,, up to p, for a range of values of J are 
a9) | tabulated in Table 2. 
: TABLE 2 
Partition factors: loads on chairs for unit load over one chair: 
parts in 10,000 
art ’ Do Ps Ds Ps Ps Sum o| remainder 
we 20 5 ° © ro) 
Ve 7 30 0 ° ° ° 
2 35 2 ° ° ° 
3 9 I ° ° 
7 15 oO ° 
207 31 5 I 
j 109 I 10 3 
123 169 51 S I 
15 130 145 79 ) 
I 123 112 ol 
> 
When the chairs are represented by an equivalent continuous elastic 
ut foundation (1) the deflexion due to unit concentrated load at x = 0 is 
y = (B/2K )e-P*(cos Ba + sin Bx), (38) 
where 8 = A/4E/, and K is the load per unit length on the foundation 
+ required to produce unit deflexion. Since a length A of the foundation is 
to represent one chair, 5 |/Kh and 8p h®/48 EI: hence 
i’) J 8./dp 12/(Bh)*. 


The load equivalent to that on the chair under the load is then the integral 
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of Ky over the interval — th to 4h, that is 


Dy = 1—e-1#* cos(4 Bh) = ¥Bh{1—jh(Bh)?+ 2(Bh)®—sbo(Bh)*+...}. (39) 


When Bh > 0, J + o and both formulae (37) and (39) yield p, > (3/47). 


10-—~- > 






0-9} 

Rail or s aced | 
0-8; deat | 
0-7} | 


06 } {Rail on continudus _ | 


elastic foundation 


(Timoshenko) 
0:5} T 1 
™N 
0-4} + + + + + + | 4 
6-= deflection of a single chair under unit load 
0:3} 6, = deflection of o single bay of rail simply 
supported at adjacent chairs under unit 
load at mid-span 
0-2 + ‘ + 
Flexible <—————— Ra | >. St iff 
O-1}— Stiff <——Chair ——=» Flexible 
“|, 2. 7 2. oe 
Values of J=6-/5p 
Fic. | 


For moderate values of J, however, the estimate of p, afforded by formula 
(39) may be about 8 per cent. less than the value given by formula (37); it 
is interesting to notice that Timoshenko recorded that ‘The difference 
between the calculated and the actual weight never exceeded 8°,’ (ref. 1, 
footnote, p. 408). The comparison between the two formulae over the 
median range of J is shown in Fig. 1. For low values of J (high values of B) 
the value of p, indicated by formula (39) oscillates about unity; this portion 
of the curve corresponding to Timoshenko’s analysis is not shown in Fig. | 
because clearly the conception of continuity of support must fail when 
either the chairs are very stiff or the rail is very flexible. 


(b) Deflexion under the load 


Another particular case is when « = 0; then formula (34) reduces to 


S _ ) | oe > 
y = Sp [ J(2+eos*)+y*(3+y y—yoos?) cos Mé dé, (40) 
7 J 2+-cos 0+ 1J(1—cos 6)? 
0 
where now y = 2v(1—v). 


The integral (40) may also be evaluated in algebraic terms; for instance, 
when v = }, so that the load is midway between two chairs, the deflexion 











YR Ul 


This | 


load 
valu 
large 
larg 
whic 
any 

redu 
tota 
chai 


(Cc) . 
A 


to 


































THE DEFLEXION OF A SPRUNG RAIL 


yp, under the load (M 0) is 


39) (YR/8x) = IPo+{h(3)!+ (14d) /[2{24- 3414 J} (41) 
J) This deflexion is thus always greater than the deflexion p, 5., when the same 






0-9 
4 
0:8 
wa 
00 Tends to i (SV 
0-5245 
0-5 
} ° Both tend 
Yc /Op to (3/4J)4 
oe 
. 6-= deflection of a single chair 
6, = deflection of a single bay 
3 of rail 
Y-= deflection of continuous rail 
. when load is over a chair 
‘ Yp= deflection of continuous rail | 
/ ~~ when load ts midway 
we between two chatrs 
| = ds fo zero 1 t t 
lula 
): 11 ; 7 r 4 
nce ~ Values of J=6-/b, 
a Fic. 2 
the load is directly over a chair. The comparison between these two principal 
fp values of the deflexion under a moving load is shown in Fig. 2; when J is 
tion large they tend, of course, to the same value (3/4./)!5,.. However, when J is 
‘4 large, the ratio of Yr— Po % tO Op, that is (Yp 87) J po: tends to 1(3 4J)}, 
her which itself tends to zero. The fluctuation of deflexion as a load passes over 
any given size of rail with any given chair spacing may thus always be 
reduced to any desired extent, but only at the expense of increasing the 
total deflexion to ten or a hundred times the deflexion of the rail between 
) 
chairs. 
(40 (c) Deflexion at chairs when load is at midspan 
\ third quite special case is when v = « }. Then formula (34) reduces 
to 
6. [f (11+-cos @){cos M6 + cos(M—1)6' ,, 
nee, | y ‘ t ate ( y 5 dé. (42) 


; 7 8(2-+-cos 6) +-2J(1—cos 6)? 
K10n bs 
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When M = Oor | this represents the deflexion at the two adjacent chairs 


_ 
when the load is applied midway between them, and similarly M = 1 or? ; - 
represents the deflexion at the next pair of chairs. Alternatively, by the a 
reciprocal principle, formula (42) also represents the deflexions midway A= 
between chairs when the load is applied over a chair. The integrated form y= 
of formula (42) for MV = 0 or 1 is } y= 
y, 1 3(K+1)?—2 . ‘= 
§. K3(K+1)?—1’ (29) M : 
where K = [3{24+341+J)8}. Sn 
When J is large, y,/5, tends to (3/4)! in conformity with the limit of Sn 
formulae (37) and (41); when J +0, K > 14+4-3-!, and (y,/6,) tends to 7 
(1/8)(10—3v3) = 0-600. The function (y,/5,) decreases continuously as J is y 
Pr 


increased and it passes through the value } when J is about 3-13. 


: x, f 
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Notation 

EI = flexural stiffness of rail. 1, 8.17 

N+1 total number of chairs. 

: ‘ : 2. C.1 

h = spacing of chairs. } 0 

l Nh total length of rail. $i. 

5, = deflexion of a single chair under unit load. L 


c 
5, = deflexion of single bay of rail under unit load midway between two 
chairs at which the rail is simply supported. 


J = 6,/8p. 

x = distance along the rail from one end. 
x = ordinal denoting chair at # = rh. 

a = value of x at which unit load is applied. 

y = downward deflexion of rail at section 2x. 

Yo: Yy = Aeflexions of rail at its ends, 2x = 0 and x = I. 

Ye = 3(Yot+M)- ) 


ob (Y%—Yo) l. 
y, = deflexion at rth chair. 
Yq, = deflexion under load at x = a. 
Yy, = deflexion at x = b. 
A,,, Am = amplitudes of Fourier components sin(nza/l) and sin(mzz/l) of 


the deflexion of the rail from the chord joining its ends, 





alrs 


THE DEFLEXION OF A SPRUNG RAIL 


p any integer, including zero. 
db ra/t. 

p= 2a l 24 7. 

A\=a/h. 

y = fractional part of A. 

Mm bih. 


€ = positive fractional part of A—p. 

M = corresponding integral part of A—y (positive or negative). 
Ss function of n and N, 

S function of n, N, and ¢. 

6= n/N. 


F = (1—cos 6)?/4(2+-cos @). 


S, = (N/7)48,. 
x, 8, y = ratios used in summation formulae (8 is also quoted in another 
sense from reference (1)). 
also y 2v(1—-v). 


(,, H, = functions of @ and v. 


+ eg. A = 24, » = 2} correspond to M = 0,€ = }, whereas A = 24, 1 = 2} correspond to 
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SUMMARY 
This paper traces the consequences of the assumption of body couple in classical 
infinitesimal elasticity and examines the mechanism of action across an interface 
in a continuous medium with special reference to a new multi-constant controversy 


in the theory of elasticity ‘ 


1. Introduction 


THE objects of this paper are (i) to show that elastic problems of body force | 


and body couple under given surface tractions can be reduced to equivalent 
problems of body force alone with a modified body force and modified 
surface shears, and (ii) to refute the arguments of Sir C. V. Raman and 
K. S. Viswanathan who have recently claimed the classical ideas of stress 
and strain to be untenable. They propose a new theory in which, in the 
case of general aeolotropy, there are forty-five elastic constants whilst in 
the isotropic solid three constants are needed instead of the usual two Lame 
constants. Since their arguments hinge on the existence of a stress tensor 


pq for which pg «4 gp, we have to consider the ways in which a non- | 


symmetric stress tensor could possibly arise. 


2. The stress system 

The stress system in a continuous medium arises from an examination 
of the nature of the action across an imagined interface in a continuous 
medium. We suppose that the action across an element of interface d8 
_ dS located at the element together 


with a couple G,,dS. Classical elasticity supposes the interfacial stress 


of normal fi can be regarded as a force R, 
couple G,, to vanish, and employs only the interfacial stress vector R,,. 

A consideration of the equations of motion of a tetrahedral volume 
element with three faces parallel to cartesian coordinate planes and the 
remaining face with its normal fi of direction cosines (/, m,n) relative t 
the coordinate axes then gives the equations of transformation as 


R,, = /R,+mR,+nR,, G,, = 1G, +mG,+nG,, (I) 


where R,, G,, etc., are the cartesian stress vectors and stress couple vectors | 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 3 (1956)] 
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respectively. We then arrive at the stresses pq and the stress couples y,,, 


from a cartesian expansion of these vectors, i.e. 
a a ao ne ‘ A . ms as oe i a 9 
R, = pr 1 py) ! pe k, G,, Ypx +YpyI +Yp: k (p, q ae iy Y, z). (2) 


the symbols with p = q denoting torsional couples on the surface 


‘pq 
elements on which they act, those with p 4 q denoting couples of flexural 
character. This completes the mechanism of the possible transference of 
action from element to element of the body. 

In addition to this interaction from neighbouring elements of contact 
type, there may be body force F and body couple Q per unit mass acting 
on the volume elements in virtue of the body’s presence in a field of 
intangible force, as, for example, gravitational body force or the body 
couple in a polarized dielectric solid in an electric field. We then consider 
the D’Alembertian equivalence of all these forces and couples with the 
mass-accelerations of the body. The vector sum aspect of this equivalence 
gives the body stress equation 

R, R, + — R,+ pF = pf. (3) 
Cx Cy : C2 
and the sum of moments aspect of the equivalence gives the body-couple 
equation 
é é 


G G,+2G,+pQ4+firR,]4 


CX Cu C2 





jAR,]+[kAR,] = 0, (4) 


where p and f(/;, f.,f,) are the density and acceleration respectively. These 
wre equivalent to three equations of type 


CXX CXLY CXZ 


- pF, pf (5) 
C2 oy Cz 
ind three of type 
OV rx CV ry OV rz 7. , j 
! cy -4+ pV, +yz—zy = 0. (6) 
Cx cy C2 


3. The stress-strain relations 
The deformation of the elastic body is described with the aid of the 


lenlane 
displacement vector D and the strain components e,,, or tensor strains €,,, 


given by equations of ty pe 


cu cw cl 
f € e Qe 
I y2 7] 
C2 c Y C2 
D. = u, D, = », D, = w. (7) 


Since 


Co »» these are symmetric by definition. The displacement of a 
volume element can then be analysed into a rigid body displacement of 


translation D and rotation w = }curlD, together with a displacement, 
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dependent upon the strains and vanishing with them, which we term the 
pure strain of the volume element. 
We now divide the stress tensor pg into symmetric and antisymmetric 


parts o,, and «,, according to the usual scheme pg = o,,+-«,,., where 
Ong = 4(pq+9p) = Fgp> tng 4(pq—qp) —= —CMqp> (8) } 
The typical Pyle equation (6) can then be written 
e 0 OY 2: 
Oey = = pQ,+ Y xz 4 Vey 5 OV ze (9) 


Ox by oz 


and in the absence of the interfacial stress couples, the antisymmetric | 
stress components «,,, are determined by the body couple. The symmetric 
components o,,, fallow from the stress-strain relations which we find by 
the usual strain-energy method. If the work done against the internal 
interfacial stresses pg and couples y,,, when the displacement D is increased | 
by 85D is 8W, then W is termed the strain energy per unit volume. Since | 
all the internal and surface forces and couples together with the reversed } 
mass accelerations form a null set of localized vectors, the total virtual | 
work done in a virtual displacement 5D by all the surface tractions R, 
and couples G, over S and by the body force F and body couple Q and 
the reversed mass accelerations — pf within S must be equal to that done | 
against the interfacial stresses and couples within S. Hence we find 


| SW dQ = I (8D.R,) dS + [ (8w.G,,) dS + | (8w.Q)p dQ 4 
Q S S a 


4+ [ (8D.F—f)p dQ. (10) 


a 


Using (1), we change the surface integrals to volume integrals. Then 
making use of equations (3) and (4) and the commutivity of the operations 
5 and @/éx, etc., we have on expanding the scalar products and identifying 
the integrands in the usual way, 


bW => maa *) - Yue a(S) + z= dw,(pg—qp) (p.9g,7 = 2. y;4) 


(11 


:, ., - ) 
This reduces still further, using the strain components e,,, of (7) and the 


relations (8), to 


pa 


SW = Zz Pag ttn + Vow 8Ar2t+Yyy 84 yy +Vee da,.+ 
‘i Yyz day.+ Yer da,.+ Yery 84,y; (12 | 


ow Cw, Cw, 9) } 
where é. = —, a,, = —*——? = --4, ee (13) | 
OX “ oy Oz 


zy’ 
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so that, if we are within the limit of perfect recovery, 
ow ow 


Tnq * Ypa = 


(14) 


Cog CA ng 


Assuming we are within the limit of linear elasticity this would lead to 
seventy-eight elastic constants in the case of complete aeolotropy. If, on 
the other hand, we assume the interfacial stress couples y,, to vanish as 
is the case in classical elasticity, the stress-strain relations are given by 
equations of type 


xa a yz = o,,.—4pQ,, zy = O,,+4pQ,, (15) 


where the symmetric tensor o,,, has the same form in terms of the strains 
as does the symmetric stress tensor pg when the body couple Q is absent, 
so that for complete aeolotropy we have no more than the usual twenty-one 
elastic constants and the usual two Lamé constants A and p for the isotropic 
elastic body, in which case the stress-strain relations are of type 


LX AO+ 2pe;;, y2 Hly2—dPQy, zy = Hey. + dpQi. (16) 
e,. = divD. 


where 3) a 


From (1) and (16) the cartesian stress vector is given in the case of isotropy 


by eD 
R, Adi- pVut pe + boli ‘ Q| (17) 
and the general stress vector by 
R, AS 2 + u{n A curl D]+3p[i A Q). (18) 
on 


Substitution into (3) leads to the fundamental differential equation for the 
displacement vector in an isotropic solid, using (16), as 


(A+-p)grad div D+pV*D-+ pF+ 4pcurl Q = of. (19) 
It appears from (18) and (19) that to solve for D and the stresses 


PY Fnq Tt “pq 
in an elastic solid under body force F and body covple Q per unit mass, 
with a given applied surface traction R,, over the boundary we need only 
solve for D corresponding to a symmetric stress system in a solid with a 
modified body force F’ and a modified surface traction R‘, over the 
boundary, where 


F’ = F+ieurlQ, Ri), = R,—}e[fiA Q]. (20) 


n nm 


Since (fh. R),) = (@.R,,), we see that the modification of the applied surface 
traction involves only the surface shears. 
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4. Velocity of propagation of waves in isotropic solids under 
body couple 
In Raman and Viswanathan’s three constant isotropic theory the bulk 
modulus cannot be evaluated from the experimental data for the velocities 
of propagation of elastic waves in a solid together with its density. It is 


} 


of interest to see what effects, if any, body couple has upon wave propa- 
gation in an isotropic elastic solid under no body force. 

If Q = grad¢ the equation has the same form as for zero body couple, 
since curl Q then vanishes. Taking the divergence of both sides of (19), 
putting F = 0 and considering small motion so that f = é*D/ét?, we have 


, ri, : 
1 


and the presence of body couple does not affect the propagation of com- 


Le 


pressional waves in the solid. Again, taking the curl of both sides of (19) we | 


find > 
tryst | Pouricurl 0 2 on gl 9 
= 5 w+ curlcurl Q = 0, cs = p/p, (22 
| C5 ¢ t?| 4 py 
which reduces to the usual wave equation for distortional waves, if 
curlQ = gradys so that V4 = 0. (23 


, 
In this case the experimental values of c,, c,, and p lead to the Lamé | 
constants and so to the bulk modulus, Poisson’s ratio, and Young’s modulus 
from the usual relations. 


5. Criticism of Raman and Viswanathan’s elastic theory 

A recent paper (1) by Sir C. V. Raman and K. 8. Viswanathan claims 
that three independent elastic constants are needed to describe the stress- 
strain relations for an isotropic elastic solid instead of the usual two con- 
stants which have held the field since the settling of the rari-constant 1 
multi-constant controversy in the earlier history of elasticity. The nature [ 
of their arguments is perhaps best illustrated by quotations (i) from (1), p.?, | 
on the torsion of a rod by end couples, where they say ‘we are forced ti 
recognize that the strains of a solid cannot, in general, be described solely | 
as elongations but may also include twists. Further the external stresses 


~~ 


applied to the body are couples. It follows that the internal stresses maj 
also be of the same nature. In other words the stresses in an elastic solid 
cannot be assumed to be exclusively in the nature of tractive forces but me’ 
also include torques’, and (ii) from (1), p. 3, ‘the three components of the 


angular momenta of any volume element will not vanish in dynamit | 

7 ai 
experiments or for heterogeneous strains involving rotations and whicl | 
accordingly involve torques. We therefore retain all the nine stres| 
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components in our formulation and do not make the usual reduction in 


their number from nine to six’. 

The word ‘torques’ which appears in these quotations can only be inter- 
preted as body couple Q or stress couples y,, or both, and we have to 
deduce from internal evidence in the paper which it is that they mean. 
It emerges that they really mean neither. For in a second paper (2), wherein 
they extend their theory to the elasticity of crystals, we have the statement 
that ‘when the elements of volume are small enough, the interactions 
between each element and its neighbours can be expressed in terms of 
tractive forces alone, it being then clearly unnecessary to introduce any- 
thing in the nature of couples or torques’, i.e. they agree with classical 
elasticity in taking y,, = 0. Again, on the subject of body couple they say 
(2), 70) that ‘there is no room for such a postulate since the analytical 
specification of the stresses in terms of the tensor components should itself 
suffice to describe the state of the solid completely’, i.e. they say Q = 0 
for reasons which we do not find convincing. 

Asa result of their rejection of each type of couple they must inevitably 
have from equations (6) the equality of the cross-shears, i.e. yz = 2y, ete. 
By not satisfying (6) they are ignoring one of the fundamental equations 
of motion and their theory is accordingly unsound on this score. Their 
choice of (Gu/cx, Cu/oy, Cu/ez, ev/ex, Ov/oy, Ov/éz, Cw/éx, Cw/ey, Cw/éez) as 
strain components’ is open to criticism as also are their comments on the 
role of ‘rotations’ in the classical theory. In (2) they say “The arguments 
on which the reduction in number from nine to six of the components of 
strain is based may be summarized by the statement that the elastic strains 
can be separated into what are called ‘‘pure strains’’ and ‘‘rotations’”’ 
and the latter can be ignored’. Again in (1) they say the rotations are 
ignored in formulating the stress-strain relationships’. It will be apparent 
from our brief description of the classical procedure that there is no question 
of ignoring the rotations in that procedure, they just do not come into the 
picture unless body couple Q or the interfacial stress couples y,, exist. As 
to their choice of strain components, the state of zero strain in the classical 
theory is the rigid body displacement of small translation (a,b,c) and of 
small rotation (W1, Wo, Ws) 

u A+ Gy Z— We Y, v b4 Ws LX — WZ, Ww C+ Wy Y— Wet, (24) 
whereas Raman and Viswanathan’s state of zero strain corresponds to a 
rigid body displacement of translation only, the small rigid body rotation 
(W),W», W3) Corresponding to a state of strain 


(0, Wg, Wy, We, UV, ~W1, — Wy, W, 0). 


This seems to us to be a misuse of the word strain, to say the least. In 
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our notation, but using Raman and Viswanathan’s three elastic constants 
d,,, 4,2, and d,,, their stress-strain relations are 
ow 


yz = d _--(d,,—d,.—2d 
y Ogg Cys + (44) —Ay2 Oa) ey 


b] 

xx = dy.5+(d,,—dyo)e,,- (25) 
2 ov 
zy = dy Oye (d,,—d9— 2d 44) De? 


On physical grounds, the condition for the equality of cross-shears which 


they should use is F ; 
. d,,—d,.—2d,, = 9, (26) 


which reduces their equations to the Lamé equations with A = d,, and 
p = d,,. Again the rigid body displacement (24), which should surely result 
in zero stresses from physical argument, leads on their theory to non- 
vanishing shears 

yz = —zy = (d,,—d,,—2d,,)w,, etc., (27) 
and the condition that these vanish is equation (26) again. 

We conclude that if there is any real evidence of the need for more elastic 
constants in macroscopic physics than the classical theory already provides, 
this would imply a possible existence of the interfacial stress couples y,, 
hitherto discounted in classical elasticity. It clearly has nothing to do 
with the existence of body couple. 
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EXTREMUM PRINCIPLES FOR SLOW VISCOUS FLOW 
\ND THE APPROXIMATE CALCULATION OF DRAG 


By R. HILL and G. POWER 


(Department of Mathematics, University of Nottingham) 
{Received 17 November 1955] 


SUMMARY 


A complementary pair of extremum principles are proved for a Newtonian viscous 


fluid in quasi-static flow. It is shown how these can be used to obtain arbitrarily 


lose approximations from above and below to the total rate of energy dissipation 


and hence to the drag on a translated or rotated body in the fluid. 


1. The elastic analogy 
Ir was probably Lord Rayleigh (1) who first drew attention to the formal 
similarity between the field equations for an elastic solid in equilibrium 
and those for a Newtonian viscous fluid in quasi-static flow. Displacement 
and rigidity modulus in the one correspond to velocity and coefficient of 
viscosity in the other. Also, since a fluid is generally treated as incom- 
pressible (infinitely great volumetric viscosity), the analogy requires that 
no change of volume is permitted in the elastic solid (Poisson’s ratio }). 
It is well known that an elastic state of equilibrium can be characterized 
by either of two extremum principles: one relating to the class of all displace- 
ment fields consistent with the boundary conditions on displacement, and 
the other to the class of all equilibrium states of stress consistent with the 
boundary conditions on stress. The two principles are complementary in 
the sense that the minimum of the functional occurring in one is equal to 
the maximum of the functional in the other. Under certain types of boun- 
dary conditions the common value of these functionals in the actual state 
is directly connected with that of an important unknown: for example, the 
torsional rigidity of a bar whose ends are given relative rigid-body rotations 
in their planes, or the penetration into a surface by a flat rigid indenter 
under a prescribed load. In such cases an approximate value of the impor- 
tant unknown can be obtained by bracketing it between upper and lower 
bounds corresponding to arbitrarily-chosen admissible functions in the 
extremum principles. In this way the error can be held within any pre- 
wssigned limits. Moreover, since the first variations of the functionals vanish 
in the actual state, a close approximation to the actual value of the func- 
tionals can usually be obtained from comparatively simple admissible 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 3 (1956)] 
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functions. This method has been exploited with considerable success, 
particularly in the determination of torsional rigidity (see, for example, 
Sokolnikoff (2)). 

However, we are not aware that the analogous method has ever been 
remarked on in the literature of viscous flow problems.+ An important 
group of these concerns the motion of solid bodies through the fluid at very 
small Reynolds numbers (such that the acceleration terms can be neglected), 
Here the particular unknown furnished by the extremum principles is the 
drag on the moving body. The elastic analogue would be the determination 
of the load needed to displace a cemented inclusion by a prescribed amount; 
this problem has been little studied since the stipulation of no slip on othe; 
boundaries is somewhat artificial in the elastic context. 

We shall give some illustration to show how the method can be used to 
estimate the drag in situations where a full analysis would either be out of 
the question or uneconomical in time and effort.{ For convenience we first 


establish the extremum principles in their hydrodynamical context. 


2. Extremum principles 

To simplify the notation the stress and strain-rate at a point are regarded 
as vectors in a nine-dimensional space, where the tensor components are 
taken as coordinates, and are denoted by o and e. The rate of dissipation 
of energy per unit volume is given by the scalar product oe, or by 2ue? in 
an incompressible Newtonian fluid of viscosity p. 

Suppose that a rigid body moves in a fluid which is completely contained 
within rigid boundaries at rest or extends to infinity in one or more direc- 
tions. Let the instantaneous motion of the body be specified by the velocity 
U of some point O of it and by the spin Q (no confusion need arise through 
using bold-face type for both three- and nine-component vectors). Let u 
be the actual velocity at a generic point in the fluid and € the rate of strain; 
u vanishes on the boundaries at rest, and at infinity, and conforms to the 
local velocity of the moving body at all points on its surface. We exclude 
from consideration unbounded two-dimensional flows produced by pure 
translation, since the only quasi-static motion possible to the fluid is a 
trivial rigid-body one. 

Let u*, e* be the velocity and rate of strain in any continuous incom- 
pressible flow which satisfies the given boundary conditions. The stress 

+ Indeed, Rayleigh’s analogy itself has received remarkably little attention ; for a recent 
application, see Hill (3). 

{ The method was originally indicated by the first author in 1951 at a Colloquium in the 


Department of Mathematics, University of Bristol. Prager (4) has stated extremum prin- 
ciples for a visco-plastic (Bingham) material which reduce to the ones used here in a limiting 


case. 








~~ 








wa 





distri 
the s 


throu 
both 
and 

Gaus 


whet 
More 
of tl 
the « 
varii 
insté 

N 
on t 
The 
field 


the 
righ 


not 


whe 
opp 
in t! 


loac 


He 


We 
to! 
tra 
wh 
tan 








CCe@&Ss 


mple 


beer 


tant 


ned 
lire: 
oclty 


ougt 


slude 
pur 


IS a 


con 


tress 





EXTREMUM PRINCIPLES FOR SLOW VISCOUS FLOW 315 





distribution corresponding to e* need not be in equilibrium. We integrate 


the self-evident inequality 


throughout tue fluid. Now the right-hand side is equal to (e*— €)o/y since 
both flows are incompressible, where o satisfies the equations of equilibrium 
and e, e* are derivable from velocities u, u*. Hence, transforming by 


Gauss’s theorem and using the boundary-conditicns, 
E* > E, (1) 


where E and #* are the respective total rates of dissipation of energy. 
Moreover, since the difference #* — E is proportional to the volume integral 
of the square of e*—e (the term associated with the original inequality), 
the extremum value of £* is actually an analytic minimum; i.e. the first 
variation vanishes, as was originally shown directly by Helmholtz (see, for 
instance, Lamb (5)). 

Now let o’ be the stress in any equilibrium state for which the tractions 
on the body are equivalent to a finite force F’ and couple I’ about O. 
The strain-rate €’ corresponding to o’ need not be derivable from a velocity 
field. We begin with the same inequality, but in the form 


oe’ —aeE > 2(c'—a)e, 


the fluid being Newtonian. On transforming the volume integral of the 
right-hand side, and noting that the boundaries at rest and at infinity do 


not contribute to the resulting surface integrals, there results 
E’—2(F’'U+T’Q) > E—2(FU+TQ), (2) 


where F’ and I” are the force and couple resultants (at O) equal and 
opposite to the tractions exerted on the body by the distribution 6’. Now 
in the actual state the rate at which work is done on the body by the applied 


loads F and T balances the rate of dissipation of energy, i.e. 
E = FU-TQ. (3) 
Hence, combining (1) and (2) and making use of (3), 
2(F°U+T'Q)— EF’ < FU+T < E*. (4) 
We have thereby bracketed the actual rate of work of the forces needed 
to move the body in the prescribed way. In particular if the body is simply 
translated (Q — 0) it isin effect the resultant force or drag that is bracketed, 
while if the motion is planar (U.Q = 0) it is the couple about the instan- 


taneous axis of rotation. 
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Some important deductions can be made from (4): 


(a) If a body with surface S, could be completely contained within a | 


body with surface S,, then the drag F, on S, is greater than or equal to the 
drag F, on S, (the velocity of translation and the external boundary of the 
fluid being the same). For the actual flow for S,, together with a pure trans- 
lational motion in the space S,—S,, is an admissible velocity field u* in 
respect of S,. Henoe 


F,U< E*=E,=F.U o FSF 


4 
This may also be proved by the lower bound principle, though a little less 
easily. 

An application of this theorem is made in § 3 (iii). 

(b) The drag F, on a given body in translational motion through a fluid 
in a stationary rigid container C, is greater than or equal to the drag F, in 
any container C, which would completely enclose C,. For the actual flow 
within C,, together with fluid at rest in the space C,—C,, is an admissible 
field u* in respect of C,. The result follows as before. 

This theorem is illustrated by the exact solution to the problem discussed 
in§ 3 (ii). Asimilar result is true, of course, for the torque on a rotated body; 
for instance, the torque needed to rotate a cylinder in the presence of a 
plane wall increases monotonically with decreasing distance from the wall 
(cf. the solution by Jeffery (6)). 

(c) The drag on a body S tends to be increased by the presence of other 
bodies S,,..., S|, which are either fixed in position or free to move without 
restraint. For the actual flow when all the bodies are present, together 
with the appropriate rigid-body motions within the spaces occupied by 
S,,..., S,, constitutes an admissible field u* in respect of S alone. Hence 
FU < E*, where F is the drag on S when alone and £* is the energy 
dissipation in the flow when all the bodies are present. Since either the 
motions of §),..., S,, or the resultant tractions on them, are zero, E* is 
equal to the work done against the drag on S itself. The stated result 
follows. 

(d) The apparent viscosity of a fluid containing a solid suspension is a 
monotonically increasing function of the volume concentration of solute. 
This may be proved directly from equation (1) along lines similar to (c), the 
velocity being now prescribed on the boundary at infinity. 

By obvious modifications in the proof of (4) similar continued inequalities 
‘an be proved when there are body forces or different boundary conditions, 
for instance when the medium is partially bounded by surfaces where the 
tractions are specified (this situation is more natural in relation to a viscous 
solid rather than a fluid). 
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EXTREMUM PRINCIPLES FOR SLOW VISCOUS FLOW 


3. Illustrations 

(i) Translated sphere in unbounded fluid. The exact solution is well known 
and was obtained by Stokes (see Lamb (5)). We use the problem as an 
illustration for the sake of clarifying the details of the method of approxi- 
mation. 

Take spherical polar coordinates relative to the centre of the sphere and 
its direction of motion. The simplest velocity field satisfying the equation 
of incompressibility, the no-slip conditions u* = U cos, v* = —U sin@ 
on r = a, and vanishing at infinity, is 

u* 2a at? v* asin@ 
0 (= “eos 7=-— 


The corresponding stream-function is 


a 


lar : 
yb* kl a*( —1)sin?0, 


Ous* . ra) * 
ia rsin@év* — — oyp 


where r?sin 6 u* ——, ~~ 
06 er 


This may be compared with the actual stream-function 


an ee 
ys = }Ua?|— ——})sin’6. 
a r 


In view of the incompressibility the rate of dissipation of energy is calculable 


as 


E = 8p | | (e2+<€,e9+e§+eF9)r? sin 6 drdé, 

a 0 
where (e,, €9, €,g) are the (r, #) components of the strain-rate tensor, namely, 
eu u le iC Cu ov. a 


€ ; €A + : € - = 
or r r 00 a rod r or 


9 


With the above choice of u*, v* it is found that Z* = S2yaU. 
For the lower bound take the equilibrium state 





, ss cos 6 ‘ p.UBa* sin @ 
Oy, pl (a . - — - 


Tr — ’ 
6 r > 
2 2r3 


r - 


with all other components zero. Here a, 8 are dimensionless constants whose 
values are to be chosen to make the bound as large as possible. To calculate 
E’ note that the corresponding strain-rate components are 

€, = 0,/dp, €9 = —o,/6p, €19 = To/2p- 
We obtain 2F’U—E’ = snpaU[6a—(a?+o08+36?)]. 


[his expression is greatest when a = 3°, 8 = —}8, and its value is ?7paU. 
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The actual drag F is therefore such that 
2 < Fi/mau < ¥. 

Thus F/zpaU = 5-97 with a maximum error of 0-25, as compared with 

the actual value 6. 

(ii) Translated cylinder in a concentric cylindrical container.+ Let a be 
the radius of the moving cylinder and 6 that of the container. Take cylindri- 
cal coordinates relative to the common axis and the direction of motion, 

For the upper bound assume a stream-function 


» «af A(b\2"  B/r\2 \? ry? 
b* l/r sin @ a ft 1 : (; +C : +D(;) ? 
2 n\r n h r b 


, 1 oyb* Oys* 
where u* p ; ve — — ; 
r © cr 


For given n the constants A, B, C, D are determined by the conditions 


u* U cos6@, v Usin@ on r= a, u* = 0= v* on r = bD; n itself 
being chosen so that the upper bound is as small as possible. If F is the 
drag per unit length, we find that 


F _ 4 n?) 
47 2n 


+ (k"1—1)(4C0+ BD)+k(k"-1—1)(AD+ BC) + (1—n2)A Blog k, 
where k = 6?/a?. The least value is obtained by letting n decrease towards 
zero, with due regard to the dependence of A, B, C, D on n. There results 


(k2"— 1)(A?+ B?)+ 3(k?—1)(C?+ D?)4 


5 < MA— Blog k-+(k—1)fMk-+ (C249) —(AD + BO)), 
4rpl : 2 
where 
A B C 
: , ot yk—2(k—1)]-. 
k+ bi , = D = [(k+1)log k—2(k—1)] 
‘ 2__727-1 
Whence I = < View b 6 a?)-* 
4zpl a 624 aa 


In the limit the A and B terms in the stream-function take the form 
r( E+ F log r)sin @. 

The above bound is actually the exact value of the drag, as may be 
recognized from the solution by Frazer (7) or, more conveniently, from 
that by Stevenson (8) to the elastic analogue of this problem. 

(iii) Translated torus in unbounded fluid. Consider the solid r = asin@ 
obtained by revolving a circle of diameter a about a tangent; the solid may 
be regarded as a closed torus. Let the body be moved in infinite fluid in 


+ This is, of course, not strictly a steady-state problem but can still be treated on the 
quasi-static basis if the motion is sufficiently slow. 
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EXTREMUM 


the direction of its axis of symmetry. The exact solution does not appear 
}be known. We obtain close bounds by the use of Theorem (a) alone. 
The body can contain a circular disk of radius a moving broadside on. 
The drag 16uaU’ on the disk (Lamb (5)) is thus a lower bound. The body 
in be contained within an oblate spheroid whose generating ellipse has a 
major axis of length 2a and a radius of curvature at the ends of the major 
xis equal to a. This requires the eccentricity to be 1/V¥2 and the minor 
ixis to be of length .2a. The drag on such a spheroid is 4v27zyaU (Lamb (5)) 
nd is an upper bound (Stokes’s formula for a sphere is a less good bound). 
Thus the drag on the torus is 5-37zpaU’ with an error not greater than 0-27 


nthe numerical coefficient. 
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THE ACCURACY OF THE METHOD OF 
CHARACTERISTICS FOR PLANE SUPERSONIC FLOW 
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[Received 31 August 1955] 


SUMMARY 
The mechanism by which errors are propagated through a computation is analysed, 
and the resultant effect of the various types of error determined. A procedure is 


outlined for planning computations to achieve specified accuracy with minimum 
labour. 


1. Introduction 

THE progressive nature of characteristics computations suggests that errors 
may have serious cumulative effects. Previous studies of related problems 
have yielded the order of magnitude of the error introduced in the individual 
steps of the computation (1), and have established the conditions for the 
convergence of the numerical procedure as well as the order of magnitude 
of the resultant effect of the errors introduced under those conditions (2). 
For a practically useful estimate of the effect, however, an understanding 
of the mechanism of error propagation and a more precise measure of error 
size are required. Because of the lack of these it has been necessary to 
expend a great deal of labour on the computations to ensure sufficient 
accuracy. The following is a study of the error problem for computations 
of two-dimensional, steady, homentropic, irrotational, supersonic flow of a 
perfect gas. While this is sufficiently specialized to admit an explicit 
solution with direct practical applications, the methods employed are not 
restricted to the treatment of this single type of computation. 

The computation procedure is briefly as follows. The characteristic 
differential equations are replaced by approximate difference equations 
which define the mean slopes of elementary Mach line segments. Then, 
starting from given initial conditions, a step by step process is carried out: 
in each step a pair of the difference equations is solved simultaneously to 
determine the point of intersection of two segments. The complete com- 
putation consists in determining a set of points of intersection, the ‘net- 
work points’ associated with a network of two families of ‘plus’ and ‘minus’ 
Mach lines. 

At each step ‘rounding errors’ (due to the finite number of digits carried) 


[Quart. Journ. Mech. amd Applied Math., Vol. IX, Pt. 3 (1956)] 
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and ‘mean slope’ errors (due to the approximate nature of the difference 
equations) are introduced. Additional errors are introduced when the 
initial conditions of the computation are known only approximately. All 
these can be represented by errors in the position of network points and 
will be considered in this form first. 

The position of the network point determined in each step depends on 
the positions of two adjacent network points determined in two previous 
steps (§2). The effect of any particular error will therefore be felt 
only in a part of the network, the ‘domain of influence’; to determine the 
effect there, the mechanism of error propagation must be studied. On the 
other hand, the resultant effect at any point will be due only to errors 
introduced in another part of the network, the ‘domain of dependence’ ; 
to determine this resultant effect, a study of the accumulation of effects 
is also required. 

A plausible conjecture concerning the nature of the domain of influence 
of an error, which reconciles the intuitive idea that ‘errors will be propa- 
gated along Mach lines’ with the wider domain indicated by the Uniqueness 
Theorem, is that the principal feature should be a direct propagation along 
the two Mach lines emanating from the point where the error is introduced, 
and that the effects in between should be comparatively negligible. This 
conjecture is confirmed in § 5. 

An important result (§5) is that the domain of dependence of a 
point cannot even approximately be identified with the two Mach lines 
converging to the point. In general, the errors introduced at network 
points between the lines make an appreciable contribution to the resultant 
effect; a contribution which is of the same order of magnitude as that from 
the errors introduced on the lines themselves. For the special case of the 
statistical effect of rounding errors, however, it is found (§7) that 
the domain of dependence may indeed be approximately identified with 
the Mach lines. 

The explicit results obtained render possible not only a check on the 
accuracy of completed computations, but also a planning in advance for 
specified accuracy with minimum labour. A computation in the design of 
a supersonic nozzle (3) has in fact been replanned to check on the usefulness 
of the results, and it was found that the accuracy desired for that compu- 
tation was attainable with about a sixth of the labour actually expended. 

This work was carried out in partial fulfilment of the requirements for 
the degree of M.Sc.(Eng.) in the University of Sydney, and the author is 
indebted to the University for a research grant. The author is grateful to 

Dr. R. E. Meyer for his guidance in the preparation of this paper, and to 
Professor A. V. Stephens for his encouragement. . 
Y 
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2. The numerical method of characteristics respect 

The derivation of the characteristic equations, and their adaptation t 
practical computations, are well known (for example, 3). For two 
dimensional, steady, homentropic, irrotational, supersonic, flow of ;! 
perfect gas, we have sugges 


} 





d 
y tan(@—,) | 








dx on any plus Mach line 
‘ to the 
6+-w x = const. 
dy 
“ = tan(@+p) | > ; 
dx on any minus Mach line 
@—w = P= 7 
where @ is the angle made by the velocity direction with the x-axis, p ji 
the Mach angle, and w is given by ) 
q 
cot pu 
w(t) f dq, 
. q 


a, 
where q and a, are the velocity magnitude and critical sonic speed respe 
tively. 

For computations the differential equations (1) are replaced by the f 
difference equations, 


3. TI 
y(a, B)—y(a, B—AB) 

t (a, Af | (% 
x(x, 8)—ax(a, B—AB) foe Oe OP) | (9 : | 
5 - step, | 

y(x, B)—y(a— Aa, B) 
; t_(a—Aa, B, A of the 
x(a, B)—a(a— Aw, B) ‘ ae ener 


where a, 8 are the characteristic coordinates of a point («,y), and Aq and —_ 
AB are the step sizes. The difference quotients represent the slopes of the (Fig. 
two straight lines (elementary Mach line segments) joining the points (a.f) | C™P 
and (a«,8—Af), and («,8) and (a—Aa,f). The quantities ¢, and f are } lines 
approximations to the mean slopes of the Mach lines connecting the had . 
respective points, whence mean slope errors are introduced. } The 


’ : ° : ‘ : (a,.B 

To determine the mean slope errors, consider a minus Mach line in an} ‘0 
: ae | are I 
ideal flow. Suppose that «, y and the characteristic parameter « are related I 


. - ont £ iS als 
along the line by x = x(a), y = y(a). Let A and B denote two adjacen 
points on this line, with coordinates 

+I 


Hy = U0), Ys = Ylay) and ay = x(aytAa), yp = y(ay+Aa), 


simila 
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| respectively. Then. for example, the mean slope approximation,t 


tion dy 
r ty f tan(? Max, bAa (") 
: aa x4+) Aa 


suggested by Temple (4), implies an error in slope, 





} 9 se ora , as , 

dy Yn—Y4 Aa®(ex\-*[@x cy Cy Cx ‘ 
At_ = | } : a it Ss 
AL} » +A Lp—Xq 24 \ea Ga® Gx Ca® Ga 
4 + , a 
to the first order. 
} 
IS, | S\a%e.s) 


espe 





3. The propagation of errors 
Consider a computation in which only a single error is introduced. Each 
“| step, excepting that in which the error is introduced, consists of a solution 
/ of the pair of exact difference equations corresponding to the approximate 
equations (2). The error is introduced in the determination of the coordi- 
vand| lates of the network point (a 9,8»), thus displacing (a, By) from O to O" 


of t Fig. 1): it is assumed that OO’ is small compared with OR, RS, etc. The 
3 (a computation is represented by the network shown by full lines. The dotted 
tare) limes show how the network from (a ,8,) onwards would have appeared 
 the| had no error been introduced; they represent the ideal form of the network. 


The figure shows that the error OO" at (a»,B,) is propagated not only to 


in al 19 Py), (%», By), ete., along the Mach line « ¥9, Where the resultant errors 

elate wre RR’, SS’, ete., respectively, and similarly to points along B = fp, but 

jacent | '8 also propagated to (a,,8,), (a1,B.), ete., where the resultant errors are 
UU’, VI et respectively. 

If the a »roximation ¢ tan $|(0+ u)4+(0+)p,)] were chosen instead, a treatment 


vould ld similar results, because ¢(@+- j2)/€a 4(1+dy/dw), ete. 
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From §2 the Mach line segments connecting adjacent points of the} The 
network are parallel to the corresponding ideal segments to a high} §<1 
order of approximation.t The error propagated from one network point| Fro 
to an adjacent one is measured by the displacement of the connecting; § = f 
segment from its ideal position, or by the intercept made by the corre 
sponding ideal and displaced segments on a segment of the other Mach 
line family. Changes in error therefore take place only at network points 
For (a, 82), for instance, the incoming error along RS is represented by 
RR’, and the outgoing errors are SS’ along ST and S’S’ along SV. Nov 
SS’ marks the continuation of the propagation along the original Mach | 
line, and is, in fact, only a modification of RR’, the change being due t; 
the difference in slopes of the segments LR and MS. Outgoing errors of 
this type will be called transmitted errors. On the other hand, S”S’ marks 
a ramification of the propagation; S”S’ is propagated along a Mach line 
intersecting the original one. Even though the magnitude of S8’S’ js 
governed by RR’ and the angle between R’S’ and MS’, the error S8’S' is 
introduced solely because the two segments MS’ and S’V’ are not collinear 





and so outgoing errors of this type will be called parallax errors. 

It has been shown that the mechanism producing transmitted and 
parallax errors depends only on the slopes of the Mach line segments, that 
is, only on the ideal network, and this yields an important lemma: Th 
propagation of any introduced, or transmitted, or parallax error is independent 


of the presence of any other error. 





| 
4. The equations of propagation 
. . . . . . ) 
The characteristic variables, « and £, associated with the points of a! Putt 
ideal network, are adopted as a curvilinear system of coordinates in thi 
flow plane. to tl 
The error at any point O (Fig. 1), is defined by two vectors, OA ani 
OB, in the directions of the pair of minus and plus, elementary Mach line | 
segments, respectively, that converge to O in the course of the computation | .. . 
4 a . . . . ° Dal 
OA and OB are written a and b respectively, and the positive directiot | es 
for both is taken to be that making an acute angle with the local stream } 
direction. Fig. 
' neit’ 
+ e.g. the slope of a straight, elementary, minus Mach line segment is, from equation (3 ine 
ni equi 
} 4 Aa? G, 
tan(@+ p) F+Aa tran 
where F is a function of the characteristic coordinates and step size alone, and G depend T 
on the position of these in the flow plane as well. For a displacement that is small compare 
with the mesh width, the change in G is AG O(Aa?). Therefore the corresponding chang? | a se 
in slope is . i mn, 
ae A[tan(@+ )] O(Aa*), : Phe 


which is small compared with the mean slope error (itself O(Aa?)). 
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The step sizes, Aa and Af (in radians), are assumed to be O(8), where 
l. 
From Fig. 1 the parallax error introduced onto the minus Mach line 


8, at (a», 8.) is given by 


S" 8’ {sin(S’SS”) sin(/ YS" S)|SS’. 


oe ) ° 
} | } 
— | 4 





10 T 








3 
Mach 
Number | 
M 2 
and a(a-f) 
f(y) _ | _d, 
Fp) min 
l(a gy? | 
160 
0 


























Ab and SS’ a, this becomes 
Ab = —o(a—f)a Aa (4) 


to the first order, where, from (5), 


Putting S”S’ 


o(a—f) [2(0+-2)/@x] cosec 2u = —(1—N)cosec 2p, 
N b(y t 1)sec*u. 

Similarly, it is found that the parallax error introduced onto a plus Mach 
line is Aa = o(a—f)bAB. (5) 
Fig. 2 shows that o(a—f) is O(1) in the practical range of (a—f) where 
neither sonic nor vacuum conditions are approached, and therefore, from 
equations (4) and (5), the ratio of a parallax error to the corresponding 
transmitted error is O(8). 

To determine the change in anerror in the course of transmission through 
a series of points along one Mach line, draw SP parallel to RR’ (Fig. 1). 
Then the transmitted error at (a 9,8.) is given by 


SS'/SP sin(S’ PS) sin( SS’ P), 
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that is, (a+Aa)/a = sin(LRS)/sin(LRS+S’'SP). 


As 6>0, (LRS) > 2u and (S’SP) + [é@(6+p)/e8] dB NdB, and the 
equation tends to 


1+da/a sin 2u/sin(2u-+ N dB) = 1—Neot 2u dB. 
Thus the change in an error in transmission along a plus Mach line to the 
first order is given by 
da/a N cot 2u dB. (6 
Similarly, for transmission along a minus Mach line, 
db/b N cot 2uda. 7 
Since Cp/ex = —Op/eB 1(1—2N), 
equations (6) and (7) can be integrated along plus and minus Mach lines 
respectively to yield 


aay I (#) f (Ho): (8 
and b by = f(H)/f(Ho). (9 
wheret f(p) {| (y—cos 2)’ sin ‘y+ Du}! Y Dsec pu}. 


In Fig. 2 f()/f(H) min IS plotted against (a—f). The ratio remains bounded 
in a computation because of the sonic and vacuum limits, but it can 


problems, which is characterized by the rapid attenuation of the effects 
of an error. 

Now consider an elementary segment of a plus Mach line. It is con- 
venient to regard it as beginning just after one network point (a, B—As). 
and terminating just after the next («,8). Suppose the error at the begin- 
ning is a(a,8—Af). To allow for other possible errors suppose that at («,8 
a parallax error, o(a—f)b Af, and a new error, Aa,, are introduced. Then, 
from (8) and the lemma of § 3, the error at the end of the segment is 

a(a,8) = | f(u)/f(u—Ap) ja(a, B—AB)+Aa,;+o(a—f)b AB, 
or a(a,B) = a(a,B—AB)+-Aa,+o(a—B)b AB, 
where a= a/f(p), Aa, = Aa;/f(u), and b=b/f(p). 


Therefore, the total change in error along a plus Mach line segment is 


Aa Aa ;-+-o(a—B)b AB. (10) 


Similarly, it is found that the change in error along a minus Mach line 


given by 
segment is given by Ab = Ab,—o(a—f)a Ac. (11) 


+ This result is the same as that obtained by Meyer (5) for the growth of first-order 
disturbances along Mach lines. 








nevertheless be large. Therefore, though an error remains of the same | 
order of magnitude in transmission, it may be markedly magnified or | 
diminished. This is in sharp contrast to error propagation in elliptic | 
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cons 


If, as 5 > 0, a;(a, B) | » Aa,| ; , and b,(a, 8) = [> Abi, ia tend 


to differentiable functions of « and f, then (10) and (11) tend tot 


ca/eB éa,/eB+-o(a—B)b, ob / Cx 0b ;/Cxa—o(a—B)a. 


u 


5, The solution of the equations of propagation 

A method of successive approximations is employed to sum (10) and 
11) simultaneously and so determine the resultant error at any point («, f). 
For convenience the field of summation is formally extended over the 
smallest Mach quadrangle containing the domain of dependence of («, 8), 
because the shape of the domain itself may differ from case to case. Thus 
the quadrangle is defined by the pair of Mach lines through («, 8), and two 
other Mach lines, « v, and 8 = By, say. Each of the successive summa- 
tions is over only one characteristic variable, and a prefix to the summation 
sion is used to denote the constant value of the other characteristic variable 
during the summation. 


For a first approximation to a, put 
x 

b(a, B’) B’ S Ab. 
— 


in (10) and sum along « const. Then 
5 B x 
A,(x, B) “> Aa > Aa; } > a( x—f’)AB B> Ab;. 
8 Xo 


Po Po 


Next, substitute a,(«’,8) for a in (11) and sum along 8 = const. to obtain 


bo(a, B) >) Ab 8 Ab BY o(a’ —B)Aa o> Aa;— 
. B x 
p> a(a’— B)A x > a( x’ — 8’) AB B’> Ab,. 
Xo Bo Xo 


Continue by substituting b,(«,8’) for b in (10) and again summing along 


. = const., and so on, to obtain finally 


Ua, B) >) Aa, 


+> o(xa—P’)AB B’> Ab 


> o(a—f’)AB Bd a(x’ —B’) Ax v> Aa; 
x Bo 


p’ x 
“> o(a—f’)AB ¥ o(a’—P’)Aa oS o(a’—P” AB p’> Ab;+ 


x Bo Xo 


n equations which are the special case, éa;/€B 
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b(a, B) = B> Ab,— 


x B 
— B> o(a’—B)Aa o> Aa;— 


Bo 
B a’ 
_ aS of B)Aa a’¥ a(x’ —B’)AB p> Ab; + 
Bo Xo 
p’ 
+ BX of a’ —B)Aa aS a(x’ BAB eS o B’)Aw «oS Aa;+ 
Bo Bo 


+.... (13) 

In (12), for example, the first term in the series gives the effect of the 
errors introduced on the plus Mach line through (a, 8), while the remaining 
terms give the effect of the remaining errors, which, but for the introduction 
of parallax errors, would not have contributed to the error at (a, 8). From 


(12), 


u < (Th/8g)A > (2-"A)"+[A/(W8,)]B > (2-*A/2)", a 
n=0 n=0 
where A = fl I, W?2, W = |o(a’ —B’) max? 
lr, = |a—apl, B= |B—Bol, 
on = |Aa min? 5g Ap min? 
A = |Aa; max? B= |Ab; max* 


Hence the series for a in equation (12) is absolutely convergent. Moreover, 
since the series in (14) converges rapidly, it is to be expected that the series 
for a in (12) should, under suitable conditions, be rapidly convergent also. 

Finally, from (12) and (14), a = O(e/5), where ¢ is the order of magnitude 
of the introduced errors. As expected, the effect, 


SA 
> aj; 
of all the errors introduced on the Mach line « = const., is O(e/5). However, 
it is found that the sum of those terms in (12) giving the effect of errors 
introduced off the Mach lines through (a, B) is also O(e/8). 
A similar consideration of the convergence of the series in (13) yields 
similar results. 
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For the special case of only a single error, given by Aa,, Ab;, and intro- 
(13) reduce to the following: 


for x Xs a Aa,, 
for a ~ dQ, = Ab; o(a—f,)AB- 
Aa; >) a( x— B’)AB o(a )—P’ )Aa : 
Ab; >) a(a- B’)AB aS of a’ — B’) Aa o(a’ = Bo) AB+ 
(15) 
for B Bo, b Ab,, 
and for 8 / a b - Aa, a( Xo B)Aa— 
Ab; >) a(a’ —B)Aa o(a’—By)AB+ 
x B 
Aa;8> o(a’—B)Aaa> o(a’—’)AB o(ag—P’)Aa+ 
(16) 


Equations (15) and (16) give the effect of an error in its domain of influence. 
At any point on the two Mach lines through the point of introduction the 
error induced is O(e) and is given directly by the transmission relation, 
8) or (9). At any point between those two Mach lines, however, the error 


induced is only O(de). 


6. Mean slope errors 
The position error at the end of a minus Mach line segment AB resulting 
At_, is 


A B cosec 2yu At_ cos?(6+-p). 


from an error in slope, 


Ab 


(17) 
To obtain the mean slope error, (3), in a form more convenient for practical 
computations, length parameters h, and hg (5) are introduced, which are 
lefined by 


CXL) Ca h ' COS( 0 


62/08 = 
ey/ep 
The length parameters, and their derivatives @h,/éa and @hg/éB, depend on 
the flow 


+L), hg cos(6- ), 


éy/ex = h,sin(@+p,), hgsin(6—p). 


pattern, and are here assumed to be bounded. From (3) and 


17), therefore, 


Ab, (Aa3/12)| p(a—B)h,—o(a—B)@h,/Ex] = (Aad/12)f, (18) 


Where o is given beneath (4), and 


p(a—B) = 4 cosec 2p €?(6-+-p)/éa? = 1N(1—2N)sec*u. 
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Similarly, the mean slope error introduced by the approximation t_ cal 
t, = tan(6—p)g.ia8 and - 
can be shown to be 
Aa, = (Ap> 12)[ p(a—B)hg +o(a—B)éhg ep | (AB3/12)d. (19 
The effect at any network point of all the mean slope errors introduced = 


during the computation is determined by substituting the values given 
by equations (18) and (19) into both (12) and (13). Since ¢ and & in 
(19) and (18) are O(1), the individual mean slope errors Aa; and Ab, are 
O(8). From §5 the resultant errors are therefore O(¢«/5) = O(8?), ie} And 
an order larger than the individual errors introduced. Even so, the resultant} stanc 
errors can be made as small as desired by reducing the step size. The 





‘ 
number of network points increases in proportion to 5-*, however, and | 
hence the labour of computation would be increased by a factor rough) , 
equal to that by which the resultant error is reduced. l 

Th 
7. Rounding errors the a 


Unlike mean slope errors, rounding errors are random, and cannot be} toa 


specified individually by the particular flow pattern and step size. Never- 
theless, their effects can be predicted statistically. 
Suppose the number of network points in the domain of dependence of | thet 
(a,8) is S, and suppose the number on the Mach lines « = const. and 
8 = const. are L and M respectively. Let rounding errors given by Aa,, A) 
be introduced at each network point. The effect at (a, 8) due to each of the By ( 
introduced errors is given by 


X = Aig, a = Aa,, 
x 5%, a = P,Aa,;+ P,Ab,. 
Bo: b = Ab,, whe 


—~ 
LO 


i 


B+ &., b Q, Aa; + Q,. Ab;, 


where, from (15) and (16), P,, P,, Q,, and Q, are functions of «, 8, and step | Sine 


size. Then, from the lemma of § 3, the effect due to all the introduced devi 
errors is L s-L ; 
a= Gart 3 BAat Rab), | lV 
t=1 f=] (20 
Me S M 
b= ¥ (Abimt S (Q1Aa,+Qz Ab,), 


and these are in fact equivalent to equations (12) and (13). 
The effect of rounding errors, A’t, and A’t_, in the mean slopes ¢, ané 








(19) | 


duce | 
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1_ can be found as follows. At any point, A’t, and A’t_ are related to Aa; 
ind Ab, respectively by 

Aa, hg AB cosec 2p A’t, cos?(6—p) ‘a 

Ab, h, Ax cosec 2 A’t_ cos?(@+-p2) TT at... 
Therefore the standard deviations of the classes to which the corresponding 
Aa, and Ab, belong are 
SD(Aa,) ‘sh oe |SD(A't,) 
SD(Ab;) )|SD(A't_) 


And from (20), since the Aa,, po are independent random variables, the 


(21) 





standard deviation of the resultant effect is given by 


(22) 


WV 


L S—L 

SD(a))? = > [SD(Aa,) + ¥ {PISD(Aa,) + PZ SD(Ab Pe | 

t=] r=1 
, S M 

SD(b)? = > [SD(Ab;)\ + S {Qi SD(Aa,) + OF SD(Ab,)F 
m=1 s=1 


The rounding errors Ax and Ay, introduced at each step by rounding off 
the a, y coordinates of the network point calculated in that step, are related 
toa pair of errors, Aa;, Ab;, by 

Aa, Ax cosec 2 sin(@—p)+ Ay cosec 2 cos(@—p), 
Ab; = Ax cosec 2 sin(é+.)— Ay cosec 2 cos(6+-p), 
+] 
—s Aa, = Aa,/f(n) v, B)Ax+ U;(«, B)Ay, 
Ab; Ab, /f() ao B)Ax+-Ve(x, B)Ay. 


By (20), therefore, we have 


L S-—L 
a= > (U,Ar+U,Ay)+ > (X, Ar +X. Ay), | 
l=1 r=1 (23) 
us S-—M ; ’ 
b= > (VAxrt+KAynt+ > (¥,Ar+¥, Ay), 
m=1 s=1 


where 


A,= BLU,+8), X,-= BU,+B),, 
y Q10; 02M, y, QU, I Q2 bp. 
Since the Av and Ay in (23) are independent random variables, the standard 


deviations of the resultant effect are given by 


I 
SD(a)P? = ¥ {UYTSD(Ax) P+ UY SD(Ay) Ph 
1 
S-—L 
¥ {X3[SD(Azx)}?-+ X3{ SD(Ay)}, 
. — . (24) 
SD(b)P = ¥ {Vi[SD(Azx)}?4 re (Ay) 304 
m=1 
Yi[SD(Ax) + YISD(Ay)P}, 
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Now for both types of rounding error, the number of independent random 
variables contributing to the resultant error, given by a and 3, is large, 
O(8-*), and each variable is distributed with a zero mean. Hence the 
resultant errors, whose standard deviations are given by (22) and (24), have 
approximately normal distributions, by the Central Limit Theorem (7, for 
example), and have zero means. 

The orders of magnitude of | SD(a)}* and [SD(6)}* in (22) are 

OfS af (Aa 1 | 7 | (Aa J lencon + [(Ad,); — , 

O{s [(A5;)in mean ay [(Aa,)3 }mean+[(A5;)3 meant» 
respectively, since P,, P,, Q,, and Q, are O(8), from $5. These expres- 
sions, and the similar pair corresponding to (24), show that, much in contrast 
to the result found for mean slope errors, reducing the step size will increase 
the effects of rounding errors! Observe also that the effect of rounding 
errors introduced off the Mach lines through (a, 8) will be negligible if those 
errors are only no greater in order of magnitude than the rounding errors 
introduced on the lines through («, B). 


Obviously, rounding error effects may be greatly reduced, with only a | 


relatively small labour penalty, by carrying additional digits. If this is 
considered together with the fact that on increasing the step size the effect 
of mean slope errors is increased by only as much as the labour is reduced, 
it can be deduced that for minimum labour the effects of rounding errors 
should be small compared with the corresponding effects of mean slope 
errors. 


8. Initial errors 

Initial errors may take many forms and general rules for calculating 
their effects cannot easily be given, but they can be satisfactorily treated 
by the methods already presented. It can be shown, however, that if the 
length in the characteristic plane of a typical Mach line in the domain of 
influence of the initial errors is small compared with unity, the effects of 
the parallax errors may be neglected in calculating the resultant error at 
any particular point, provided one at least of the initial errors lies on a 
Mach line through that point, and provided the particular Aa; or Ab; trans- 
mitted along that line is not much smaller than each of the other Aa; or Ab; 
Under these conditions (and in fact whenever the parallax errors can be 
neglected), the effects of initial errors are independent of the choice of step 
size for the computation. 


9. A note on planning computations 
Although h,, hg, @h,/@x, and @hg/@8 cannot be found exactly without 
computing the flow itself first, satisfactory approximations to them call 
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usually be obtained from the prescribed values of h, and hg on the 
lines from which the computation starts.t| Such approximations are 
readily found when h, and hg are continuous and slowly varying, but 
even when this is not so (5) rough estimates can be made. Any discon- 
tinuities in h, and hg can be accounted for, since they persist along Mach 
lines and therefore must appear on the initial lines. The formulae (18) and 
19) should only be applied to segments on either side of such discontinuities. 


mations must be used (5). Every branch line intersects an initial line, 
however, and it is therefore known in advance where such modifications 
will be needed. 

A procedure for planning a computation so that the desired degree of 
accuracy is obtained with a minimum of labour can be summarized as - 
follows. Firstly, estimate the effects of the initial errors on the assumption 
that parallax errors can be neglected. At least, this indicates whether the 
desired accuracy is attainable.t Next, find the smallest number of steps 
for which the resultant errors due to initial and mean slope errors combined 
lie within the desired limits. For this, estimate the values of ¢ and ys (equa- 
tions (19) and (18)) over the network; obtain a tentative scheme of step 
sizes by inspection; and, adopting this scheme, check on the resultant 
errors: this indicates any modifications necessary. Lastly, determine by 
inspection the minimum number of digits that can be carried subject to 
the condition that the effects of the rounding errors are to be small compared 


with the correspondir:g effects of the mean slope errors. 


For the nozzle of (3), for instance, it can be shown that it is sufficient to assume h, 
nd ch and hz and @h,/é8, to be constant along plus and minus Mach lines respectively, 
and to be given by their values on the initial lines. 

t If the effect of these parallax errors is indeed negligible (§ 8), this preliminary result 
will define a new tolerance for the remaining errors, and the effect of the initial errors 
need not be considered any further. 
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SUMMARY 

Small vibrations of a conservative homogeneous system, doubly coupled statically 
are analysed. The system considered consists of a ‘et of flywheels with spring 
connexions between them and with a spring between each flywheel and the ground, 

The characteristics of the coefficients of the polynomials are considered and the 
recurrence formulae are given ; and it is shown that when these formulae are applied 
it is not necessary to develop the Lagrangian determinants. 

The periods can be determined by means of binomial equations or by a method 
of finite-difference equations of the second order. Their magnitudes can easily be 
determined graphically. 

Mathematical expressions of sums of combinations of certain trigonometrical 
functions of argument 7 are also obtained. 


1. Introduction 
THE kinetic and potential energies of a holonomic system which performs 
small vibrations about its configuration of stable equilibrium are homo- 
geneous quadratic forms in the generalized velocities (7) and coordinates (q) 
respectively, with constant coefficients (1), namely, 
T = \4A{q}, -V = 4aCiq}. (!) 
Here A = (a;,) represents the square inertia matrix, C = (c;,.) the stiffness 
matrix, {q}, {q} the column matrices, and q, q the row matrices (2). 
By Lagrange’s equations, the system of differential equations of motion 
may be written in the matrix form 
A{G}-+C{q} = {0}. (2 
To solve this system we assume that {g} = {r}(exp twt), where {r} is the ampli- 
tude matrix and w the circular frequency of vibration. The problem is 
then reduced to one of characteristic values, and the period equation is 
A,,(A) C—AA| = 0, (3) 


where A w” is the characteristic number. 


2. Frequency equations of some particular systems 
When the system has no inertia couplings, and the elastic terms are 
derived from two sources, i.e. when a; = 0 for i ~ k, (2) becomes 
D{q}+Ciq} = {9}. (4) 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 3 (1956)] 
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where D = (a,,) is the diagonal inertia matrix and C = C,+C,. By a 
INC) suitable transformation this problem can be reduced to a special problem 
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D In the special case of a homogeneous system, in which case the coefficients 
of the matrices are given by a;; @, Cy = 6, Cr, c’, (3) becomes 
A (A) AI—P—Q 0, (5) 
art where P = D-!C, and Q D~'C, are the dynamic matrices. 


ror the three characteristic cases which correspond to linear oscillations 


of several coupled masses (Fig. 1), or to the torsional vibrations of light 
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shafts with several flywheels (Fig. 2), after reducing the stiffness of the 
spring to the equivalent torsional rigidity, the corresponding matrices P 
and Q are found to be 


p—p 09. 90 fH) 2p —p 0. 0 ¢O 
—p 2p—p . 0 0 —p 2p-—p . 0 0 
Po-| . . ... .|, PO= 
0 0 0 2p p 0 0 0 2p p 
0 0 0 —p Pp 0 0 0 .—p 2p 
(6 
2p —p 0 0 q 
—p 2p —p 0 q 
Po = : Q 
if) 0 O . 2p—p , 
0 0 0 .—p p q 


(p = cla, q = c’/a). 


The determinantal equations AW? = 0, A® = 0, A = 0 corresponding 


to these three cases can be found by the use of formulae of regression, 
obtained by simple expansion of the determinants in terms of principal 
minors; moreover, A‘) = 0 and A“) = 0 areexpressible in terms of A) = 0. 
Thus ™ , “ me z : 

An” = (A—2p—q)A,,” :—p?Ay 2 = 9, Ay” I, \/ 
and the determinants yield 


A@ = (A—g)A, = 0, 


n—1 
(c) — (b 2A (bd) == 
An = (A—p q)Ay” :—? AG i 0. (5 


It follows that the case in which both ends are clamped (Fig. 1 (b)) is 
the basic one, since for the other two cases (7) may be deduced from the 
equation for this case. 

The polynomials and their coefficients are given in Table I. 

Since the coefficients B,_, are completely determined it is not necessary 
to develop the determinants (7) and (8). 

For example 


A®) == 1; A® = A—2p—q = 0, 
AS) = (A—2p—q)AY— p? = A®—(4p+ 2q)A+ (3p?+-4pq+-q?) = 0, 
A®) = (A—2p—q)AY)— p?A””’ 


= M—(6p+3q)A?+ (1Op?+ 12pq+ 3q?)A— (4p? + 1L0p?q+ 6pqg?+q3) = 0. 
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For n = 3 the coefficients are: 


r= 0, B, = 1, 


2+-s 


1. (6+8 
r=1, B= (-1) Z 3)( s |p = —(6p+ 39), 
s=0 
(=z, B, = (—1)? (ae as = 10p?+ 12pq- 39’, 


3 P 
r=3, B,=(-—1) 2 (75) = —(4p*+ 10p*q+ 6pq? +9"). 


~ 


From (8) one obtains 
A®) _ (A—q) A?” 
= M—(4p+ 3q)A?+ (3p?+ 8pq-+ 3q?)A—(3p?q+ 4pq?+43) = 0. 


The coefficients are: 


2 +\ /2+. 
r= 1 B, = — Ss (; )( ‘\pow = —(4p-+ 39), 
. 3 —S§ s 
2. /3-+8\/1-+s8 
r=2, 8,= > (; )( ‘pas = (3p*+-8pq+3q”), 
7s” 
2 1316 
r=3, B=— z iy 2) = —(3p*q+4pq?+-¢°). 
= | 


In the same way we find 
AY) = (A—p—q)AP—p?A® = d2?—(3p4 2q)A+(p?+3pq+q?) = 0. 
For the more special case when c’ = c, i.e. when p = q, the coefficients 
B,,_, of the polynomials form a diagonal series of numbers whose differences 
n-r : g 
are Ar = (—1)"(3p)’, Ar+i — 0, (9 


where r is the index of the row, and may be computed from the formula 
B a (1 A)? B® = BO 4 J fii ABO 4 _(e—-r Ar B® (10 
= +4) i. ie a l ae ae , 0° 


The determinants (7) and (8) are now of the form: 
A) — (A—3p)A 1— pA”, = Y, 
A@ = A—p)A® , = 0, ( 


A® = (A—2p)A ,—p?A®, = 0, 
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and the coefficients B,_, of the polynomials f(A) = 0 are connected by 


P 
relations (Bo? J. = [Bo ” n—pl Bo) ie 
) | an vie = [BY e+ 2 Bel, 1° 


Their values are given in Table IT. 


(12) 











' 
4 TABLE II 
n|B,| B B, B; si B, BR | | B, 
q°) , l 1 | Pp 
2 ] —4p | 3p? 
} 3 1 7p 14p? 8p 
4 | -10p 34p? —46p* 21 p4 
5 l 13p 63p? —141p* 145p* | 55p® 
6 | l6p 101 p? 320p3 534p* | —444p% 144p° 
7 l 19p 148p? 610p% 1431 p* 1905p5 1331 p® 377p’ 
V, 8 | 1 22p | 204p? LO38p* | 3160p4 5878p5 6512p 3926p? 987p® 
l l -3p 
2 l —6p | 8p? 
3 l — 9p 25p? — 21p* 
4 l —12p 5lp? | —90p* 55p* 
5 l Lip 86p? | —234p3 300p* 144p5 
6 18p 130p? 480p3 951p* | —954p> | 377p* 
7 21p 183p2 |—855p? 2305p* |—3573p> | 2939p* |—987p? 
8 | 1 24p | 245p? 1386p* | 4740p* |—10008p> | 12707p* |—8850p7? | 2584p8 
) l 1 — 2p 
} 2 l —5p 5p? 
3 l —8p 19p? | —13p% 
! 4 l Llp 42p* 65p* 34p* 
j l4p 74p? 183p3 215p* 89p5 
6 | 177 115p? —394p3 717p* | —654p5 233p* 
7 l 20p 165p? 725p 1825p4 2622p 1985p® 610p’ 
8 23p | 224p? 1203p* | 3885p* |—7703p5 9134p* |—591lp? | 1597p* 








- 0), 

¥ cients) °° Determination of characteristic numbers 

erences| Substituting € = A—2p—q = [(p?/n)+-7] in (7) we obtain 
A, nA, 1 (p? n)(A,, 1 A, 2): 


nula | Ce Varies from 1 to n—1, the sequence of equations becomes 
TT ld 


4,=1,A,—7 (p?/n), Ag—nA, = (p?/n)?,..., A,-1—7A,,-2 = (p?/n)"> 


and we obtain 


A, nA 1 (p* n)’ or 7 "A. —7'-"A,_, = (p 7)” — g2n. 





Summing the n equations of this form obtained by substituting n = 1, 2,...,n 


} We have the geometrical series 
A nr g*—. } 


] Ya’ =a— ' 
ia v=1 a—1l 
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From A, = 0 one obtains iis 
, a = i, 
° ) 7 vat —_ vit Vir 
1.€. é= r( 1!) — plexp +-onp-— 2p cos 
n p n+ 1 n+ n+ 1 
and the characteristic number is 
vit OF 
A, = q+2p{1-+ cos = q+4pcos* ; (13 
n+] 2(n+-1) 
The characteristic values belong to the interval q < A < 4p+q. They 
are given in Table III for all three cases; here A, > A, > >A, 
TaBLeE III 
Case a b c 
v l 21 
A, | A A q+2p\l+cos—7 q+2p\1+cos T q+2p\1+ cos - 
|v > Aveal @ 2p| cos ~ n) 1+ 2p\ COs - ;7) 1+ 2p\ ey, 
‘ ae L a ae n t i ee n 
AAS <A q-+ 2p(1 cos : 7) q-+ 2p(1 cos x) q 2p(1 cos 7 7| 
v v v+l n n+] 
1 a ae n—1 n ay ee n v By eissces n 


The method of finite-difference equations of second order may be used 


here also (3). The differential equation of torsional vibrations of the it! 
disk (Fig. 2) is given by 


J, ¢, [Cy-3(9,—9p-—-1) —Cx( 9x41 — 9.) |—ey, 9, (14 


Assuming that the motion is harmonic with amplitude A, and with « 


circular frequency w, we substitute 0,. 


Cp_-y Appt (Cp_1 +0, +0, —AI,) A, 


A,,cos wt and obtain 


-C, Ay, = 9, (15 


, 


or, for the homogeneous system (J, = J, c, = ¢, ¢, = ¢’), 
2 2 = 
A,-4—(2+u—v?)A,+Azy, = 9, (It 
where u=c/e, v? = Jw?/e = JAIc. 





Substituting in (16) A, = 


cosh z 1, z is an imaginary number, z = if, where i = ,/( 
have coshi8 = cosf, i.e. 


2cosB = 2+u 


v= = 2+-(q/p)—(A/p); 


J, and the characteristic number is 
A = q+2p(1—cosf). 
The coefficient 8 must satisfy the boundary conditions. Putting 


k 0, 


here p=c/J,q=c’ 


lLn,n-+l 


exp(kz), we obtain 2coshz = 2+u—v*. Ii 
) 
1), and we 
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I 
n (15) we obtain 
(Cy—Ady) Ag—Cy A, = 9, (18) 
Cy Ag+ (eg+c+ce’—AJ)A,—cA, = 0, (19) 
; cA,,_,+(e+¢,-4 vif AJ)A, Cr A aad (), (20) 
’ Cy A n T (C,, NJ, 4)A n+l 0. (21) 
(13) From (18) and (19), or (20) and (21), we obtain 
Th 1+u—v?+(n,/c)|A,—A, = 9, (22) 
; A, _,—[l+u—v?+(7,/c)]A, = 9, (23) 
- where ? “0 1 en ‘ (24) 
Wh 10 l (Co AJ) In | (c,, AJ,,) 
Substituting (22) and (23) into (19) and (20), the boundary conditions for 
> ~ mr} | 0 and f n l are 


} 
Ae e—1)A, = ¢; Anat(™ 14, 0, (25) 


Cc 
" 


Substituting for the A,’s given by A, = CcoskB+ Dsin kf, we obtain 


the following algebraical equations : 


) ° 7 
be us C(1+d¢cosf)+ DdsinB = 0; d e_}; 
the C9 
, 7] ° ° Q Qn 
(lbeosnB + ecos(n+1)8 D\ b sin nB+sin(n-4 1)B} 0: b= ] 
Cy, 
’ . 
(26) 
with ae ; 
So long as (¢ 0, D () this system of linear equations has non-zero 


solutions if the determinant of the system vanishes. This condition yields 


\ the Treque ney equation 


sin(n + 1)B+(d+)sin nB+ db sin(n— 1)B. (27) 
Let us consider the case (b). i.e. 
Jn J ; ©. ” Co c c. db us 0: 
j 0 " n" 
u then (27 becomes 
and : 
sin(n—1)p (), pb, vir/(m+- 1) 
| 
| and the characteristic values are 
| 
vir ' 
A q+ 2p{ 1—cos ; (28) 
n+] 


We can obtain in the same way the characteristic values for the other two 


cases. The latent roots are also given in Table III: but here 


A 





n° 
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From (28) it is seen that the characteristic numbers in the case (b) cay 
be easily determined graphically bydividing the half-circle of radius R = 2, 
into n+1 parts, as is shown in Fig. 3. 


f-- iets jth Chott 








Fic. 3. 


4. Sums of the combinations 

Taking into consideration the relations between the roots (Table III 
and the coefficients of the polynomials (Table I1)—the generalized Vici 
conditions—one obtains, for instance, in the case (b): 


> Cr, = (—1)B,_, = (—1) S x + -(7 7 (7 Vy tg! 


r 
—— f—és 


s=0 


In the special case when p = q, it follows from (28) that 


a 9 | = a 

> 0(3—2cosw) = S ("" l—(r ”)(" () ). (29 
— y—-=§ ] 8s 
s=0 


Trigonometric formulae relating the sums of the combinations of rt! 
class with n members are given in Table IV where w has special values. 
When the homogeneous system is coupled simply, i.e. when q = 9, t! 
latent roots and the sums of the combinations are given in Table V. 
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From the above formulae for combinations we obtain the following | 


recurrence formula 


n 


— 


¥ C"™(1—cosw) = p 1)e 1—P\ SS On cosw 
— r A ( ) r p p 


p=0 


from which the sums ¥ C" cos w can be deduced. 
By means of these sums the various recurrence formulae for preceding 
sums (Table IV) can also be obtained (Table V). 


Comparing the above results one can obtain the following relations fo 


the special binomial coefficients, too: 


Case (a) 


S | n- 4 YY ( 1p 8 %(” l °) (mm 
LZ, \1+2s p 


m 


n 


J 


i 
s=0 


— (°" (r »\(" -(r ‘) 


t 
SN ( 1) 3° 1—2p 
— 
p=0 





Case (b) 


v (" ] “*) S ( 1p 3" =" " (am 
Is y p 


/ mnt 
s=0 p=0 


r—s Ss p 
s8=0 ° p=0 


The index ¢ has the same value as m only for n 
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SUMMARY 


lhree theorems concerning the values of the roots of polynomial equations are 
lerived by manipulation of the Routh—Hurwitz criteria. When applied to the 
cteristic equation of a linear system, they vield information on the dynamic 
nstants of the dominant mode. The first theorem gives a test for locating the 
ist negative of the real parts of the roots; the second gives a formula for this real 
part when its magnitude is small; and the third gives new expressions for the 
ginary parts of the corresponding root-pair. Thus the degree of stability of the 
system may be estimated, and when this is small, the formulae give the natural 
juency and the damping of the dominant mode in terms of the coefficients of the 
racteristic equation; in this condition, it is also shown how to determine the 
lynamie constants of lesser modes in certain cases. 


{An incidental discussion is also given of the effects on stability of zero coefficients 


enaract ristic equati n 


1. Introduction 

THERE has been a tendency recently to devote attention almost exclusively 
to the frequency response approach to the analysis and synthesis of control 
systems because of the power and convenience of the graphical techniques 
involved. Much design work is still done, nevertheless, by using stability 
criteria determined from the coefficients of the characteristic equation of 
the system. In these circumstances it is clearly advantageous to extract 
the greatest possible amount of information which can be used in design, 
irom manipulations of the coefficients, provided that the processes entailed 
ire not so involved as to be too cumbersome to use. 

The roots of the characteristic equation determine the nature of the 
dynamic response of the system, which is usually dominated by those roots 
with the least negative real parts, consisting typically of a pair of conjugate 
complex values. The design problem, which is to find the damping and 
lrequency of the dominant mode, thus reduces to finding the real and 
imaginary parts of this root-pair. We shall investigate this problem, using 
is a starting-point the well-known Routh—Hurwitz stability criteria, and 


shall show that considerable developments of them are in fact both possible 


ind entirely practical. 





(Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 3 (1956)] 
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2. Criteria for stability 
We first summarize some results, mostly known, which will be required 
subsequently. The characteristic equation of a linear system may hp 
written 2 
I (p) = a+, p+a, p*+...+a, p” 0, (1 
where dp, },..., a, are known real constants, and p is a complex variable 
Let the roots of equation (1) be P,, P,,..., P, and, further, let their real and 
imaginary parts be «,, as,..., x, and jB,, jPo,.--, JB, respectively. The neces 
sary and sufficient conditions that the system be stable are that 
be <= Uy Og < Bons a <= U. (2 
Alternatively, expressed in terms of the coefficients as the Routh (1, 2) 
Hurwitz (3, 4, 5) criteria, the same conditions are satisfied if 


ay 1 4, 3 } 

H, =a -0, A, 0, H. ia; Gala 0 
n—1 2 3 n n—-2 n-4 
a a 9 | 

n n-2 1 0 a a 
n-1 n-3 
Ay 1 ay 3 a, > ‘| 
a, a, 9 a, 4 - | 
QO An-1 An-3 | 
H, | Q a, a, 2 | 0 
0 0 An-4 
0 7) a, 





(3 
where the sign of a,, is to be chosen as positive.t 
A necessary but not sufficient set of conditions for stability is (3) 
y 0, a, ces t, > %, (4) 


Criteria (4) may be sharpened slightly by taking into account the effects 
of zero coefficients. The results, which are derived in the Appendix, are as 
follows: 

If any of the coefficients ap, @,..., @,_, are zero, and the others positive 
(with a, > 0), the system is divergently unstable, except in the following 
special cases, when the system may be on the border of stability: 

+ Since the equation 

An t+Gn-1 PtAn-2 Pp? +... +a) p” = 0 
has roots 1/P,, 1/P,, 1/P3, ete., and since the sign of the real part of P,. is the same as that 
of 1/F,, it follows that the criteria can also be expressed with the suffixes in (3) all deducted 
from n. That is, 
a, > 0, . ‘ > 0, etc., 


which are well-known forms. 
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DAMPING 
Case (i) When the only zero coefficients are all the first few: a, a, @9,..., @,. 
Case (ii) When the only zero coefficients are every alternate member: 
Ay; Ay, Ay,.--,4,-,, if nis odd; 


and Ba, Bas Ges. @ if n is even. 


n-1? 
Case (iii) When the only zero coefficients are all the first few and every 


ilternate member of the remainder: 


When a system becomes just unstable and in an oscillatory state, a 
complex pair of roots, which we shall take to be P, and P,, moves to the 


imaginary axis. Conditions (2) still hold except that 
4 ug = 0. (5) 
\lso, conditions (3) still hold except that (2) 
H i... = ©. (6) 


nt a” 
When a system becomes just unstable, but not in an oscillatory state, 
one real root, which we shall take to be P,, moves to the imaginary axis. 


Conditions (2) still hold except that 
x, = 0, (7) 


and conditions (3) still hold except that (2) 


3. Criteria for degree of stability 

One measure of the degree of stability of a system is the distance of the 
nearest root from the imaginary axis, i.e. from the stability border. If 
this distance is less than some arbitrarily specified amount h|, we may 
say that the system is insufficiently stable, and vice versa. Thus the 
conditions for an adequate degree of stability, when expressed analogously 
to inequalities (2), are 

te By tty Bees ~<a, (9) 

where h 0 

A method for interpreting conditions (9) in terms of the coefficients 
ly, @y,..., @, Was given in a Russian paper (6) by Cypkin and Bromberg, 
and rediscovered by Grossman (7) and Bothwell (8). Apparently the basic 
principle was first used by Frazer and Duncan (9) in connexion with the 


numerical solution of polynomial equations. Briefly, one substitutes 


p q+h (10) 
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in the characteristic equation, and so obtains a new polynomial equation 
in q 


F(q) = f(q+h) 0. (11) } 
If we write F(q) = by) +6,q+6.9?+...+6, q” (12) 
the coefficients by, b,,..., b, may be obtained by comparison with the right- 
hand side of equation (11) expanded in a Taylor series. Thus we find 
: ad } tad } "(r) } "(n) 
by = fih), b, =! =. jt at... 6, Fo. 
: 3 oo r! nN: 


For conditions (9) to hold, it follows from equation (10) that the roots 
of F(q) = 0 must have all their real parts negative. 

These conditions can be written down immediately, by applying the 
criteria (3) to equation (12); we thus obtain a series of inequalities 

=e 4 > &.. eS (14) 

where J,, J,,..., /, are determinants which are identical with those of 
inequalities (3), except that 6 is substituted for a throughout, the 6’s being 
given by equations (13). 


By way of illustration, consider the third order characteristic equation 
S(p) = 49+, p+a, p?+a,p* = 0, witha, > 0. (15) 


The criteria for adequate degree of stability are 





" LA") f(b) 
I J ) . 0, I. i 0 
=: fh) f'(h) 
3! 1! 
J hd f(r) . . (16) 
and te hy) Fy) 0 
: 3! 1! 
It is frequently unnecessary to evaluate the determinants ay Maas I, as 


the following new theorem may lead to a quicker solution. 
THEOREM I. For the equation 


Sf(p) = ag+a,p+...ta, p" = 0, witha, > 0, 


to have all the real parts of its roots less than h, it is necessary (though not suffi- 
cient) that 


f(h) > 0, fh) > 0, F(R) > Ov... f° Mn) > 0. (17) 
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To prove this, note that for f(p) 0 to be in the stated condition, it is 
necessary that b, 0. b, 0 b 0 (18) 
by analogy with inequalities (4). Result (17) then follows directly from 
equations (13) and (18). (The latter equations also yield f("(h) > 0, but 
this is a redundant condition, as it is automatically satisfied by the premise 
that a 0.) 

Theorem I may be extended slightly to cover the cases where some of 
the derivatives are exactly zero. The results, which follow directly from 
the analogy with the stability conditions given in § 2 (criteria (4) et seq.), 
ire as follows 

If any of the derivatives f(h), f’(A),..., f"-(h) are zero, the others being 
positive, one or more of the roots have their real parts less negative than h; 
except in the following special cases when the least negative real parts may, 
instead, be equal to h: 


Case (i) When the only zero derivatives are all the first few: 


f(h), f(A), f° (A). SOA). 
Case (ii) When the only zero derivatives are every alternate member: 
fth), f(a), fh)... fU-Mh), if n is odd; 
and fh), £°(R), OA), f- MA), ifn is even. 


Case (iii) When the only zero derivatives are all the first few and every 
ilternate member of the remainder: 

f(h), f(a), f(a)... LOMA), LOA), SOTA), fOMA),...; f@-MA). 
The use of Theorem I will be made clear by the following examples. 
Example 1. To find whether all the real parts of the roots of 

f(p) p? + 15p4*+ 93p3 + 303p?+ 532p+ 408 i) (19) 
are more negative than (—4). 


The successive derivatives of equation (19) may be written down at 


sight as : ‘ - 9 . - 
f'(p) 5p* 60p? +-279p? + 606p + 532 
f’(p) 20 p? + 180p? 4-558p+- 606 
(20) 
(9) 60p* } 360p 1 558 
f{'(p) 120p+-360 
Starting to substitute p 4 in equations (19) and (20), and beginning 
with the last of these for simplicity, we immediately find f(—4) — 120. 


Having arrived at a negative derivative we need proceed no farther since 
it follows from Theorem I that the real parts of the roots cannot all be 


more negative than (— 4). 
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Example 2. To find whether all the real parts of the roots of equation 
(19) are more negative than (—3). 

Substituting for p, as in Example 1, we obtain 

for—3)=0, f"(—3)=18, f"(—3) = 12. 

Here we have found a zero derivative in a sequence which is not one of 
those given above in cases (i), (ii), or (iii). Hence, one or more roots have 
their real parts less negative than (—3). 





Example 3. To find whether all the real parts of the roots of equation 
(19) are more negative than (—1). 

On substituting for p, it will be found that the derivatives are all positive, 
Theorem I is then not applicable, and it becomes necessary to use criteria 
(14). J, 4, Lg, 4, 1, turn out to be 10, 327, 19,681, 2,946,050, 294,605,000 
respectively. Since these values are all positive, the real parts of the roots 
must all be more negative than (—1). 

Actually the roots of equation (19) are —2+2j, —3, and —4+j, in 
accord with the above results. 


4. Nearly unstable systems 
For the purposes of the following theory, let us suppose h to be specially 
chosen just equal to the least negative of the real parts of the roots. Taking 
the relevant roots to be a conjugate complex pair, P, and P,, we then have 
xy = A = h. (21) 

(Compare equation (5).) Conditions (14) will still be valid, except that 
(6, 10), by analogy with equation (6), 

I I 


n n-1 


From equations (13), (21), and (22) we have 


f@ D(a) f@ 3)( x1) 


0. (22) 


(n—1)! (n—3)! 
F ie. x) | ina 2)/ x4) 








n! n—2)! | 
i. = ™ ) — 0), (23) 
0 f(a) 
(n—1)! 


| 

e eo af 

Incidentally equation (23) is the same as one derived by Orlando (4) for 
another purpose. It yields an expression for «, in the form of a polynomial 
equation. As such, it is scarcely more useful for general computation than 
the original characteristic equation. However, when |a,| is small, con- 
siderable simplification is possible. For under this condition, the general 
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formula, obtained by successive differentiation of the characteristic 


equation, 


nfo , (r+1)! n! —— 24 
j (x) ra, 1! Sr41 To Ty Om *y ’ (2 ) 
= f(x, ) = 
simplifies to 1 a,t+(r+1)a,,. 4. (25) 
7 
Substituting from equation (25) in equation (23), we find 
Any tna, Oo a,-3+(n—2)a, 2% 
}a,+(n+1).0.a, G@,-2+(m 1)a,-1% 
0 26 
it 0 Any +Na, % (26) 


the subsequent determinants we shall now restrict n to a 
4: but the methods are quite general. Equation (26) then 


For clarity in 
specific value, 


bec ymes 





9 | 
| Ag+ 404%, G@,+20,0, | 
~ ‘ | o”7 
I, Q,+5.0.0, Ag+3ag a, Agta, % | = 0 (27) 
| 
| 9 7 | 
| : Agt+404a, G,+20_ a, | 


Thus J, is a third-degree polynomial in a,, and may be written as a Taylor 


277277 37 23 

=| ,% |; | P “iF 3| — 0. (28) 
F a2 a1] 5.3 

CX, Ja1=0 2! | @a4 x1=0 3. | Cay J a.-0 


The coefficients in this series may be evaluated as follows: 


series 
[ Xy 
I(x) [J3] 0, 07 i! 





From equations (27) and (3) we find 


ag ay 
[Is]a-0 = | G: | = Ay. (29) 
ag, ay | 


Also, in differentiating equation (27) we use the rule (11) that each row is 
to be differentiated separately, and the resulting determinants added,t 


and thus 


td, “dy a3 ay a | 

( I, } 
On A, A, A ».0 3a, ay Q, Ay | (30) 

1Ja,=0 
a, a, a, ay, 4a, 2d, 
K,+ 4,+- Kg, (31) 
+ This procedure is algebraically simpler than the equally valid differentiation by columns, 
since the resulting determinants are reducible, as in equation (36). 
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where K, is obtained from H, by replacing each coefficient a, in the first 
row by (7r+-1)a,,,; A, is obtained by similarly replacing the second row 
and K, is obtained by similarly replacing the third row. 

In equation (28) we may neglect the terms in aj and a3; and substituting 
from equations (29) and (31), we then find 


x - —. 3 
1 ; (: 
K,+ A,+ ky 
We infer that the general formula for the nth order characteristic 
equation is given by 


THEOREM II. The decay-constant a, of the dominant mode of a nearly 
oscillatory system, i.e. the least negative of the real parts of the roots of its 
characteristic equation, 

Ag +A, P-+..-+ A, p" 0 
is given approximately by 


ae (33 
z ee 


replacing each coefficient a, in the j-th row by (r+-1)a,,,. 
The results of this theorem for the first few values of n are as follows. 
' a 
For n = 2: 4 = ——, (34 
2a, 
thus agreeing with the known elementary result for a second-order system. 


For n 3: 


Ay Ay 
a, 4a, A, Ay—Ap As (35) 
My 3 ae ~3\° 
3d, ay My Ap 2(a, 3+-a3) 
») 
ag ay as 
For n = 4: 
ag ay 
ay Ap A 
> & @ on 
=, As As ay Ag ay 
» € » 
2} Gy Ag Ay 3a, Ay 2|Q, GA, Ay 
» 
a, a, a, a, 2a, as 
£92 
__ 4, 4,43—Aya3—ajz ay, (37) 


» 2 
2a3(a, a.+a5— 4a, a4) 








where K; is the determinant obtained from the Hurwitz determinant H,,_, by | 
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5, Calculation of natural frequency 
When the system is just unstable, equations (5) being satisfied, the 
natural frequency of the undamped mode is 
B 
. 1 « 
i ae (38) 
“7 
Frazer and Duncan (2, 9) have shown how f, may then be evaluated in 
terms Of dy, @,,...,@,. One substitutes p = P, = jf, in equation (1), obtain- 
ng for the real and imaginary parts two polynomial equations in (j8,)* 


Ay +a,(jB,)?-4 a,(jP,)* ae 0 | 


fees of : (39) 
ind a, +45(J8,)?+45(jB,)*+.-. = 0 J 


Elimination of the higher order powers of (j8,)? using Cramer’s rule then 
ields (j8,)* as a ratio of two determinants. We shall now give a method 
for obtaining the latter determinants directly from the Hurwitz deter- 
minant H,,_,. 

The eliminant (4) (or Sylvester’s resultant) of equations (39) is a deter- 
minant identical with H,,_,. Now according to a result from the theory 
fdeterminants (12), when the eliminant is zero, the signed primary minors, 


N,, Nj,..., V,,-,, associated with any row are given by 
y T.N y ; 2(n—2)-(4 2(n—3)- *(4 2. 
N, Ng: Ngiee.Ny 12 GB)?" (By )2"- s ..  (GB JB (40) 
If we write M,, Mg...., M,_,, as the unsigned minors corresponding to 
Y,, Vg,..., V,-;, equation (40) implies the following 


THeoreM Ill. The natural pulsatance B, of a just unstable system, that is, 


j)xthe purely imaginary roots of its characteristic equation, is given by 
2 r 
e (41) 


here M, is the minor obtained from the Hurwitz determinant H,,_, by deleting 
iny row and any (except the last) r-th column, and M,,, is the similar minor 


btained by deleting the same row and the next column. 


Theorem III shows that there are many equally valid expressions for f,; 
vhich of these is the most useful will depend on circumstances. However, 
the following two solutions are of special interest. Firstly, if we choose, in 
pplying the theorem, the last row and the last two columns (r = n—2) 
ve arrive at Cremer’s (13) expression: 


» aH, 
’ H,. 


3 


(for n > 3, H, 1). (42) 


Aa 








354 A. T. FULLER AND R. H. MACMILLAN | 
Secondly, we may choose the last row and the 2nd and 3rd from lag | 
columns (r = n—3). We find a, is a common factor of both minors, and | 
thus 


2 a : (for n > 3), (43) | 


where H',_, is the determinant obtained from H, 3 by replacing each 


coefficient a, in the last column by a;_,. Note that equation (43) is simpler 


than equation (42), in that it contains fewer cf the coefficients do, a),..., a,, 
Values of 8? for the first few values of n are as follows: 


Ao 





For n = 2: Bp? — rb (44) 
9 
; ; a a 7 
For n = 3: BE xe — = -!. (45) 
a, ay 
. > a 
For n = 4: = —. (46 
ds 


. ‘ 2  @,4,—AyA, A, 4,—Aga; , 
For n = 5: Rao 22 1" (47) 
@,4,—Aj4; M3A,—A,a, 
Although in this section we have assumed the system to be on the 
border of stability, the expressions for 8, will be approximately valid for 
1 : 
a slightly damped system, i.e. for one in which |a,| is small, since the roots 
of the characteristic equation are continuous functions of the coefficients 


7 . } 
(14, 15). For a more heavily damped system, where |«,! can no longer be 


regarded as a small quantity, the expressions for «, and £, obtained from 
Theorems II and III may still be used as first approximations to their true 
values. These could then be found, in particular numerical cases, by using 
one of the various methods of successive approximation reported by Frazer 
and Duncan (9). 


6. Damping-factor of nearly unstable systems 
The damping-factor p, associated with P, and P, is defined as the ratio 
\x,/B,|, and is a measure of the resonance of the system. An alternative 
measure very frequently used is the damping ratio, ¢,, of the dominant 
mode (that is, the ratio of the actual damping to critical damping). The 

relationship between p and € is 
(2 = — p = or p* =— 2 ne 
1+ p? 1—? 
When jq,| is small, we obtain an expression for p, from equations (33) and 
any one of (41)—(43). If we choose equation (43) from the latter we arrive at 

, 
H. w (A, -3 i, 3) 


oe n= \n—3/ Ena) (49) 
Pi K,+K,+...+Ky-, 


(48) 
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Values of small p, for the first few values of n are as follows: 
a 


he 2. 1 5 50 
For ” Py 2,/(dy4,) (50) 


A,Aa,5—A,a, (/|a. 
. ‘ 1% Ag M3 3 F 
For 7 a: ps = . ; (51) 
2(4,43-+-A5) AN Ay 


Making use of the approximation H, = 0, we may put equation (51) in the 
following symmetrical form, given previously by one author (16): 


a a 


59 
Py — 5 ae (52) 
2,4/{(@, 4g +45)(Ap4_+-44)} 
; A, Ay ,—A, az —a? a, - 
For n 4: py — : (53) 





2(a, dg +a3—4aya,),/(a, 43) 


7. Equations for the lesser modes 

It is interesting to note that, when the system is on the border of stability, 
the characteristic equation can be factorized and the factors may be 
expressed very simply in terms of do, 4)...., a,. This result follows from 
Trudi’s theory of elimination (12, 17). The characteristic equation may 
be written in the factorized form as 


f(p) = C.(p?+B}).g(p) = 9, (54) 


where C' is an unimportant constant, f? is given by equation (43), and g(p) 
is the determinant obtained from the Hurwitz determinant H,_, by substi- 
tuting for the Ist, 2nd.,..., (n—1)th elements of any one column the terms 
p)"-?, (—p)"-4,..., (—p), 1. The relative importance of the lesser modes 
of a nearly unstable system can be assessed by study of the equation 
y(p) = 0, which might be called the subsidiary equation of the system. 
A further interesting result (4) is that the Hurwitz determinants asso- 


ciated with this subsidiary equation simplify to H,, Mg,..., H,,_», i.e. to the 


n 
first (n—2) Hurwitz determinants of the original characteristic equation. 
It follows from this and equations (6) and (8) that if r roots lie on the 
imaginary axis, the r determinants H, to H,_,,, are zero.T 

Various alternative expressions for the quotient g(p). for the first few 
values of n, are: 

: ( Pp) Xo Ay ees 

For n 3: g(p) oc p- : (55) 

a, ay 


This result has been given previously for r even (18). 





A. T. FULLER AND R. H. MACMILLAN 
lag (—p)? .| 


Forn = 4: g(p) ale (—p) 9 


i's l a, | 
a aA, a. 
- 91 3 03 ; . | 2 ~p 
oc p*-+— p+ —— = = p®?+20, 8. p+ B3, (56 
a, a,a, 


us 
rm 
| 


1 /fa,a. a 
where 2=s . 3) and pf? ——*.. (57 
2A) \do dy a, BF 


a a, (—py 


‘ a; a3 (—p)? 7 


For n = 5: g(p) (58 
- G@ (—p) ay] 
| 
a; I a, | 
oc ps4" p* 4 Oy (43 My- a, 45) 4 Ay (43 U4 — Ay 45) 
a A;(A, @4—Ay) 4s) ;(A, Ay—Ay 45) 
a a 
g®4.— pt 41 9 2 (59 


a; a; B? a; B? 


8. Systems having large damping 

The results given in the present paper for zero-damped systems are com- 
plementary to, and in some respects analogous to, results previously :lerived 
for systems in which the dominant mode is critically damped (16, 19, 20), 
(For the former systems, a, = a, = 0 and B, = B, + 0; and for the latter, 
xy x» ~ 0 and B, = f, = 0.) In support of this analogy, note that for 
systems which are completely aperiodic (8, = B, = ... = B,, = 0), criteria 
have been given (21, 22) which are similar to the Routh—Hurwitz in- 
equalities. 

We may also mention at this point that various other criteria are avail 
able, e.g. for finding whether the roots are restricted to a sector of the 
p-plane (10). Marden (15) gives an extensive bibliography. 


9. Conclusions 

Theorem I provides a necessary but not sufficient condition which is 
useful in finding whether all the real parts of the roots of the characteristic 
equation are more negative than a specified amount. Theorem IT provides 
an approximate formula for the least negative of these real parts, valid for 
nearly unstable systems. Theorem III, and in particular equation (43), 
provide new formulae for the natural frequency of nearly unstable systems. 
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These results, which involve simple manipulations of the Hurwitz deter- 
ninants, may be adapted to yield expressions for the damping-factor. 
similar methods give the characteristic equation in a simple factored form. 
Some of the expressions relating to just-stable systems are analogous to 
those for just-aperiodic (or quasi-critically damped) systems. 


APPENDIX 
The effect of zero coefficients in the characteristic equation 


Consider the characteristic equation 
a) +a,p+a,p*?+...ta,p" = 0 (60) 


4s having some of its coefficients zero. We assume that a, + 0 since we wish to deal 
with an nth order equation, and that the non-zero coefficients are all positive, so 
nstability cannot be immediately inferred by criteria (4). 

If the zero coefficients are changed to small negative values, the system is then 
lefinitely unstable from criteria (4). But the roots of any polynomial equation are 
known to be continuous functions of the coefficients (14, 15). Hence, if we now 
hange the coefficients back to their original zero values, the system must either 
remain definitely unstable, or move on to the border of stability. Consider the latter 
ilternative. The roots of equation (60) must then include (a) zero roots, or (6) purely 
maginary roots, or (c) both zero roots and purely imaginary roots. These cases 
may be treated separately as follows: 

a) If there are (r+ 1) zero roots, equation (60) may be written in factored form as 


prtt(C,+C, p+C,p?+...+C,_- 1p") = 0. (61) 
Since the remaining roots have their real parts negative, application of criteria (4) to 


e polynomial in the brackets of equation (61) shows that C), C;,..., C,_,_, are all 


positive. Hence, comparison of coefficients in equations (60) and (61) yields 


Sua > & Gis > Ge & > O (62) 
and further « omparison gives 
& = @, = 4,=.. = 6, = 0. (63) 
b) If there are k pairs of conjugate imaginary roots, equation (60) may be written 
n factored form as 
DP \(yo+ p?).-(ye+ p? (dy +d, p+d, p?+...+d,_o,p"-*) 0, (64) 
where Vy 0, y2 6.5 > ©, (65) 


Multiplying out the left-hand factors in equation (64) we obtain 
(D)+ D,p?+ D,p*+...+ Dy, p™*)(dg+d, p+d, p?+...+d,_,p"-**) = 0, (66) 
where, from expression (65), 
A>), B> 6... Be> & (67) 


Also, since the remaining roots have their real parts negative, application of criteria 
4) to the polynomial in the right-hand brackets in equation (66) shows that 


d 0, d Bhs Rie > SD (68) 
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When the brackets in equation (66) are multiplied out, it follows from expressions 
(67) and (68) that there are no zero coefficients in the resulting polynomial, except 
when there is only one term in the right-hand bracket, i.e. when k 
even. In the latter circumstances, comparing coefficients in equations (60) and (66), 
we obtain 

a, =a, a; we = Gy, = 0, (69 
and, using expressions (67) 


a, > U, a, > 0, a, > 9...., 4, > 0. (70 


(c) If there are (r+ 1) zero roots and k pairs of conjugate imaginary roots, it follows 
as in case (a) that 

Ay a, —a,=>..=—a,= 0. (71 

On dividing the characteristic equation by the factor p’*", we are left with an equation 

similar to the type treated in case (b). Hence there are no further zero coefficients 


except when 


— _ - _— ew) 
Arie = Bis = Anse +5 = Ap = Y, (72 


and then Gus > BU, Guz > 0, @ b Diy ey > D. (73 


r+5 


Summarizing the results of the three cases, we find that the distribution of the 
zero coefficients can only be given by expressions (62) and (63); or by expressions 
(69) and (70); or by expressions (71), (72), and (73). 

Conversely, when a characteristic equation has some zero coefficients (the others 
being positive), the system may be on the border of stability, if the zero coefficients 
are only (i) all the first few, or only (ii) every alternate member, or only (iii) all th 
first few and every alternate member of the remainder.t Apart from those rather 
special circumstances the system is divergently unstable when there are zero 
coefficients. 


+ Note that cases (i), (ii), and (iii) do not quite correspond respectively to cases (a), (b), 
and (c). 
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PROPAGATION IN A GYRATIONAL MEDIUM 
By LL. G. CHAMBERS (University College of North Wales) 
[Received 9 June 1955] 


SUMMARY 


The electromagnetic properties of media such that D = eE+€H, B = (E+ uH 
are discussed. It is shown that each field component obeys the equation 


VF +iw(C—£)V x F- w? (we —CE)F - 0 
and that such media are accordingly doubly refracting. Integral forms for Maxwell's 
equations are also discussed. 
1. Introduction 


TELLEGEN (1) has suggested the possibilities of media whose constitutive 
relations are of the form 


3 

D, = 2 (alt $a), (1.1la) 
3 

B, = 2 (Se E,+ Ha Hh), (1.1b) 


a time factor exp{iwt} being assumed. Very little is known concerning the 
properties of these media when the quantities €,,, f,, do not vanish. (When 
they do vanish, the constitutive equations are those of ferrites the properties 
of which have been discussed fairly fully in recent years (2).) For the sake 
of simplicity, it will be assumed that the media are homogeneous and 
isotropic. The constitutive relations are thus given by 


D = «E+ €H, (1.2a) 
B = CE+-uH. (1.2b) 


A particuiar case of such a system is an optically active medium (3) whose 
constitutive relations are 


D = cE-+iwgH, (1.3) 

B= —iogh+-pll, (1.3) 

where g is the polarization coefficient. Solving (1.2) for E and H, one 
deduces that (we—LE)E = pD—EB, (14a) 
(we—CE)H = «<B—CD. (1.4b) 


If pe—lé = 0, the medium is ideally gyrational and this case will be dealt 
with in another paper. An attempt will be made here to provide the field 
theory for gyrational media which is already available for normal media. 
Let it be noted in passing that p, «, €, ¢ may all be functions of w. 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 3 (1956)] 
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2. Description of field by potentials 


Maxwell’s equations are 


VxH—iwD = 0, (2.1 a) 
Vx E+iwB = 0. (2.1 b) 
It follows immediately that 
V.D= 90, v7.B—0 (2.2) 
and from (1.4) that V.E 0, V.H = 0. (2.3) 


Using these relations and feeding equations (1.3) into (2.1), it may be shown 
that each of the four field vectors B, D, E, H obeys the vector differential 


eel VF + iw(f—2)V x F+-w*%(we—CE)F = 0. (2.4) 
Writing 7 = iw({—€&), k* = w*(we—E), this becomes 
V°F+7V x F+/°F = 0. (2.5) 
Now (2.1 b) is satisfied identically by 
B= ¥xA, (2.6a) 
E -VV —iwA. (2.6 b) 
Substituting in (1.2 b), 
H WV <A+EUV +iwfA}, (2.6) 


and substituting in equation (1.2a), 


vxas for 4 tina. (2.64) 
B be B 


Substituting equations (2.6c) and (2.6d) in equation (2.1a), it may be 


us 
us 


D eVV iweA r 


verified that 


V2A iw(C £)V <A w (je —CE)A ViV.A+i(pe—LE)V} ams (i, 


(2.7) 
On application of the divergence operator to equation (2.6d) we find that 
0=V.D = (®- (nr: iwW.A). (2.8) 

be 


It follows, by comparing equations (2.7) and (2.8), that if the vector and 


scalar potentials are related by 


V.A+iw(eun—E&l)V = 0 (2.9) 
the vector potential obeys the same differential equation as the field 
ee VA+nVxXA+MA = 0 (2.10) 


and the scalar potential obeys the equation 
VV+KhV = 0. (2.11) 
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Utilizing the relation (2.7), the fields may be expressed entirely in term; 
of a potential A; e.g. 


E = —WV—iwA = —22[V(V.A)+4A]. (2.12 

Clearly a similar set of relations may be obtained by commencing with 
D = Vx A*, (2.134 

H = VV*+iwA*, (2.13b 


which satisfies (2.1 a) identically. 

By proceeding in a similar manner to the above, it may be shown that 
V*, A* satisfy exactly the same differential equations and relation as 
V, A and equation (2.13b) may be rewritten 


ua ‘e(V(V.A*)+KA*] (2.14) 


It will be noted that if » = 0, the problem is formally the same as that for 
a normal medium (i.e. one of non-gyrational type). 


3. Condition for zero loss 


The energy dissipated within a surface S is given by (4) 


re [ 4(Ex A).dS, (3.1) 
s 


where the tilde indicates that the complex conjugate value is to be taken. 
Now , 7 ; 
(Ex H).dS = | V.(E < Al) dr 
S Vv 
— | (A.(Vx E)—E.(VxA)] dr 
V 
= [ iw(E. D—B.A) dr 
a 
= iw | [E.(¢£+2A)—A.. (wH+E)] dr 
vr 
= iw { [2\E\*—p/H|2+ (€—QE. Al dr. (3.2) 
Vv 
Now for the medium to be lossless, » and @ must be real, and & = , that 
is, w and ¢ real and € = &. When these conditions are satisfied 


and k? = w®(we—LE) = w*(we—EE) 


are both real. 
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4, Plane waves in unbounded medium 

Let us consider solutions of Maxwell’s equations which depend only 
upon distance measured in a particular direction (which for convenience 
will be taken as the z-axis). Then 


Va k— (4.1) 
O74 
and Maxwell’s equations become 
kx © tall = 0, (4.2 a) 
C2 
k ». nl iwD =— 0. (4.2 b) 


C2 
Each of the field vectors satisfies 
V*F+7V x F+/F = 0, 
oF oF 


or using (4. l RB a ae nk > = + k?F = (), (4.3) 
oz" C2 
Now 0 k. | k > =| = —iwk.B. 
Oz | 
Thus k.3 = 0 (4.4a) 
and similarly k.D = 0, (4.4b) 
and consequently k.E = 0, k.H = 0. (4.4, d) 


Thus all the field components are perpendicular to the direction of propa- 
gation and the field may be termed Transverse Electromagnetic. 

Suppose now that the fields vary as expf{iyz}, ie. there are harmonic 
waves moving in the negative z-direction with phase velocity c = w/y. 


Maxwell’s equations become 
k > iyE 1iwB = 0, 
k x tyH—iwD = 0, 


or cB = —kxE, (4.5 a) 
cD = kxH. (4.5 b) 

It follows that B.E = 0. (4.6.4) 
D.H = 0. (4.6 b) 


These two orthogonality relationships take the place of but coalesce into 
the one relationship (5) E.H = 0 when the cross coefficients € and { vanish. 
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5. The doubly refracting property 








Cer 
Let F represent any one of the field vectors. As F is orthogonal to the | follow 
direction of propagation, it must be of the form very § 
, 2) 1: . - 5 9 } 
F = (F,i+F,j)exp{tyz} (5.1) | ) 
, , a : : , from 
if the wave is propagating in the negative z-direction. homes 
rol 
Equation (4.3) becomes ie 
ones This | 
—y"(FLi+F,j)+h(£,i+F,j)+7.iyk x (F,i+F,j) = 0, since 
the exponential factor having been dropped. may 
Resolving, we have Us 
F(k?—y?) — inyF, == @, (5.2 a) 
Finyt+F,(k@—y*) = 0. (5.2b) 
Eliminating F, and F,, we get 
| e—y? iny| _ 
| iny = | Thes 
or (k?—y?)?— y*y? = 0, (5.3) 
n 9 “ao 2 oe 224 2 2__4}-4) 
whence ZY = —— Tey ; (5.4) 
4 
Thus, if » ~ 0, y* has two possible values and two waves can be propagated 
with two different phase velocities. The value of y? which tends asymptoti- ‘| 
. 2 , : e 
cally to 7? as n? > o is termed y*_ and the root which tends to 0 as 7? > « 
is termed y*.. Thus 2.2 : 
. nr = &, 
and hence y,y- = +k. (5.5) Suy 
r a — — F = 
lhe positive sign corresponds to two waves travelling in the same direction, 
and the negative sign to waves travelling in opposite directions. Now Th 
yk xD = ykx [e¢E+¢H] ] 
= —ewB+w&D, (5.6 a) 
ykx B= ykx[¢E+ pH] 
= —wlB-+opD, (5.6 b) 
these following by virtue of the constitutive equations and equations 
(4.6). Multiplying equations (5.6) vectorially by k and substituting back 
the values of k x D, kx B one obtains tM 
an 
(y?—w"* e+ w*é?)D+ (w2el—wef)B = 0, (5.7 a) di 
(—w*nl+- w*é)D+ (y?—w* e+ w?l?)B = 0. (5.7 b) be 
Eliminating B and D again yields equation (5.3) for the propagation con- a 
stant. 
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Certain points of interest arise in the preceding analysis. Firstly, it 
follows from equations (5.2) that the ratio F,/F, can be real only under 
very special circumstances. Secondly, if 7 = 0, then y? = k? and equations 
(5.2) become nugatory, each term vanishing identically. Thirdly, it follows 
from equations (5.7) that B and D are in the same direction, and hence 
from the constitutive relations that E and H also lie in the same direction. 
This is apparently inconsistent with equations (4.6), but is not actually so, 
since the components of the vectors are complex, and their scalar product 
may vanish. 

Using equations (5.2) possible fields are given by 


FF Fr, i Y ’ j\expfiy, 2}, (5.8 a) 
| iny, “J 


F F, be y-—" 


: expiiy_ 2}. (5.8 b) 
(  imy- *J 


F. e. F U¥+—Y-) slexpfiy, zh, (5.9 a) 
\ n } 
F_ = F414 (+7) jlexptiy_ 2}, (5.9b) 
n | 
the total value of D being given by 
F=F,+F_. (5.10) 


Suppose that F, F. f,, which is equivalent to saying that at z = 0, 


F = Fi. Let 
ytAy, y-=y—Ay. 


Then 


F = Fiiexpliy,, z}[exp{iAyz}+exp{—iAyz}]+ 


F,j exptiy,, 3] : .2y [exp —iAyz}—exp{iAyz}] 
1) 


» A 
2K, expiiy,, 2! icos Ayz-4 i727 sin Aye}. (5.11) 


7 


This will represent any disturbance moving in the negative z-direction as 
any plane may be chosen as z = 0 and the i-direction may be taken as the 
direction of F there. It follows from (5.11) that any such disturbance can 
be regarded as a wave moving with phase velocity c,, = (w/y,,) subject to 
a modulation of periodic length 27/Ay and that the polarization is elliptic. 
Suppose now that » is small. (This corresponds to the case of optically 
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active media, the treatment given above being partly a generalization of [t follows 
the theory of such media (3).) Then 





2Ay = y+ 2: = 4 a ¥: = (4h°9?-+ 9)! or 
” ” Myst+y-) nye t+y-) 
and this expression—which is exact—tends to unity as 7 > 0. If this js 
> Cas . ; . Ke These eq 
the case F = 2/,expity,, z}[icos nz+j sin nz], (5.12) 


and the polarization is circular. 
It may be noted that the same considerations apply to the other field | which 15 








components. longitud 
the TE] 
6. Wave with vanishing longitudinal component Subst 
The vector potential and the field components all satisfy 
V-F+7VxF+F = 0 (6.1) 
if F is expressed in terms of the rectangular cartesian system. Suppose | the solu 
now that, as before, the field varies as exp{iyz}. Let | 
c : ) It may 
= +> - = = = +») 
laa "5 Votiyk, (0.2) | virtue « 
F = F,+£k. (6.3) 7 Ger 
Using these relations (6.1) reduces to It is 
* 7 9 6 » f 
V6 Fy + 1k x (tyF,—V, F}+ (k?—y*)F, = 0, (6.44) 
V2 F.+7k.(V, x F,)- (k?—y*)F. = 0, (6.4b) subject 
. . . Pp ) 
Suppose now that the longitudinal component F. vanishes. Then - 
equations (6.4) become . 
V2 Fy +k x {iyF,}+(e@—)F, = 0, (6.5a) Substi 
nk.(V,xF,) = 0. (6.5 b) 
From the latter it follows that 
a : whenc 
OF, Cc F, 0. (6.6) 
cy Ox 
which is satisfied by i= x, Fk, = Fa (6.7) Clear} 
Now V.F = 0 for all the field vectors (see equations (2.2), (2.3)) and 
substituting (6.7) into this, it follows that Any 
a2 ‘.. case : 
Vex = at x = QV, (6.8) 
ox oy” Co 
and hence Ver, = 0. (6.9) subje 
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It follows therefore from equation (6.5 a) that 


iynk x F,+(k2?—+)F, = 0 (6.10) 
al (k?—y*)F.—inyF, = 9, (6.11 a) 
inyF + (k? —y*)F, == @, (6.11 b) 


These equations are consistent only if 


(k?—-y?)?— n*y*? = 0, (6.12) 


which is the same as equation (5.3), so that any wave in which the 
longitudinal component is absent, will propagate with the velocities of 
the TEM wave. 

Substituting the values for F given by equation (6.7) in equation (6.11 a) 


(k? y2) X- iny <x = (), (6.13) 
Ca cy 
the solution of which is 
x = f (inye+|h—y* Jy). (6.14) 


It may be verified that this satisfies both equations (6.11 b) and (6.8) by 


virtue of equation (6.12). 


7, General solution in powers of 7 


[t is possible to obtain a power series in » which is a formal solution of 


Vx(VxF)+7Vx F+/°F = 0 (7.1) 

subject to the condition V.2 = 0, (7.2) 

Let F = ¥7’F,,. (7.3) 
n=0 


Substituting 


x 


> o"'Vx(VXF,)+ ¥ 7" H4VxF,+ >} ky"F, = 0, 
n=0 n=0 n=0 


whence, equating the coefficients of the powers of 7 to zero, 


Vx(VxF,)+FF, = 9, (7.40) 
Vx(VxF,)+F,+VxF,_, = 9, a > 1. (7.40) 
Clearly the conditions of (7.2) and (7.3) are equivalent to 
V.F (), n > 0. (7.5) 
Any solution of (7.40) subject to (7.5) represents a possible field for the 
case » = 0, i.e. for a normal material. 
Consider now Vx(VxF,)—/"F, vx FP...» (7.6) 


subject to V.F == §. 


n 
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. It wil\ Nowais 
whence ( 


The solution of this proceeds on the lines outlined by Stratton (6) 
be seen that, if in Stratton’s equation (11) one puts 


E=F,, J=0, J*=-—F,,,, - | 
\ 


the solution of (7.6) by Stratton’s equation (19) is seen to be v 
1 The sing 
ry. r , =. 
F(t) = —- [ [F,a(r')x 0'6] ar’, (7.7) 4 and 00 
where = |r—r’|-'exp{—czk|r—r' |}. (7.8 —|u 
. . . . . y 
The negative sign has been taken because in this paper a time facto; 
exp(iwt) is assumed and the surface integral vanishes because the volume 
integral is over the whole of space. In this way, it is theoretically possible 
to build up a field which reduces to a known field for the case of » = 0. } Now 
8. Integration of field equations and on 
It may be shown, by extending the work of Stratton and Chu (6), that it } » = #. 





is possible to express Maxwell’s equations in integral form for gyrational The | 
media also. By the analogue of Green’s vector theorem 
| {(pa.[V x (V x E]—E.[V x (Vx da)]} dr ; 
. So 
v : 


= | {Ex[V(d¢a)]|—¢a x (V x E)}.dS, (8.1) 


S 





where ¢ = exp{—tkr}/r and r is measured from the element at (x,y,z) to Takin; 
the field point (x’, y’,z’); a is an arbitrary fixed unit vector. The small | 
sphere S,, 7 = 6 will later be excised. Now 


. 
Vx(VxE) = 7Vx E+E nth 
limit, 
and V x (V x da) = ak?¢d4+-V(a.V¢4). 
The left-hand side of equation (8.1) becomes 


| (ga. [E+ nV x E]—Efak’s+V(a.V)p]} dr 





) 
} . “ 
= | {a.dn(V x E)—V.[(a.Vd)E]} dr in ex 
v 4a 
-a.| [ ian dr — f V4(E.dS)|, 
Le Z 
since V.E = 0. . 

The right-hand side of equation (8.1) becomes Exal 
[ {a.[Vdx (ExdS)]+¢[a x iwB]}.dS ar 
= not. 

=a a.! ( Vd x (E x dS)+iwBd > ds). mn 
Fs | in te 
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it wi Now ais an arbitrary vector, and so may be rermoved from the equation, 














whence (8.1) becomes 


[ iwByd dr — | (Wd)E.dS = | [Vbx (ExdS)+iwBd xdS]. (8.2) 


I 
The singularity at r = 0 may be excised by a sphere S, of small radius 3, 
and sO 
| iwBynd dr i{Vd(E.dS) + Vd x (E x dS)- iwBd x dS} 
i S 


| {Vd(E.dS)+V¢ x (ExdS)+iwBdxdS}. (8.3) 


Now Vd (ik -  exp{—ikrif, (8.4) 
rj) r 


and on S, the normal is directed outwards, that is, centrewards so that 
a= f. 
The right-hand side of equation (8.3) becomes 


(i T | “exp ikd}f(E.n) dSo4 





S 


+ (ik + | “exp ikS\t x (Ex n) dS 4 iwB(5Jexp(—iks) x n ds, 


oO 
Taking the limit as 6 > 0, and remembering that 
r(E.n)+fx(Exn) n.(E.n)+nx(Exn) = E, 
the contribution of the spherical surface S, becomes 47E(z’, y’,z’) in the 
limit, and so 
47E(2’, y’, 2’) 


| iwBnd dr + [ {Vd(E.dS)+V4d x (ExdS)+iwBdxdS}. (8.5) 


j § 
In exactly the same way, it may be proved that 


4nH (2’, x’, 2") 


| iwDnd dr + | fioDdxdS — Vd(H.dS)—V¢ x (Hx dS)}. 

. S (8.6) 
Examination of equations (8.5) and (8.6) shows that there is a fundamental 
difference in the nature of the fields according as to whether 7 vanishes or 
not. If» is zero the first term on the right-hand sides of equations (8.5) and 
(8.6) vanishes, and the field within the volume V can be expressed entirely 
in terms of the fields upon the bounding surface. If, however, 7 is not zero, 

5092.35 Bb 
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this is not the case and the fields depend not only upon the boundary values 
but also on the behaviour within the volume. 

Utilizing the constitutive relations (1.2) equations (8.5) can be re. 
expressed in terms of B and D to give 


4nB(z’,y’,z') = ( iwy(uD—CB)d dr + [ iw(uD—CB)dxdS 


I Ss 
| {Vd(B.dS)+V4dx(BxdS)}, (87 
4nD(z',y’, 2’) = ( iwy(ED—eB)d dr 4 ( iwd(ED—eB) x dS 
V S 
— | {Vd(D.dS)+Vdx(DxdS)}. (84 


Alternatively, the equations could be expressed entirely in terms of E 
and H. 

The volume integrations present certain difficulties. However, by ex. 
pansion into a complete orthogonal set over V, or a Fourier transform 
it is possible that these difficulties may be overcome. 
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THE STURM-LIOUVILLE EQUATION 
EXPANSION IN TRIGONOMETRIC SERIES 
(8.7 ! By 8. C. R. DENNIS (The Queen’s University of Belfast) 
[Received 21 October 1954. Revise received 6 October 1955] 
SUMMARY 
(5.8 Ina recent paper (1) Dennis and Poots have given a particular numerical solution 
F ofa problem in heat transfer in which it is necessary to determine the discrete set 
of E ff eigenvalues and associated eigenfunctions of an ordinary differential equation of 4 
| Sturm-Liouville type. In the present paper an attempt is made to standardize ; 
} the analysis to deal with a more general Sturm—Liouville system, and to point out 


Y €X- | the possibility of a feasible general numerical method of approach. The eigen- 
‘orn functions are expressed as infinite series of trigonometrical functions, orthogonal 

wer a finite interval of integration, and the Sturm—Liouville equation is reduced 
to an infinite set of linear simultaneous algebraic equations for the coefficients of 
the series. These equations contain the characteristic parameter, and may be solved 
by standard iterative or relaxation methods. Two numerical examples are treated 
n detail. 


1. Introduction 


ly this series of papers we shall consider the solution of the Sturm—Liouville 


- 


differential system 


7 [ p(x)y’ | + (q(x) +-Ar(x) |y — 0, (1) 
y'(0) = hy y(0) ) 
y'(l) = hoy(l) J’ 


defined throughout the finite interval (0,1), where h, and h, are prescribed 
constants and A is a parameter. This system arises in many physical 
applications and, in particular, in the solution of partial differential 
equations by the method of separation of variables. It is not our purpose 
to discuss the theoretical aspects of the problem, except that we shall 
assume the well-known result that if p(x), p’(x), g(x), and r(x) are con- 
tinuous in the interval (0,/), then the conditions (2) may be satisfied only 
if A has one of an infinite set of eigenvalues A,, (n = 0, 1, 2, 3,...). Corre- 
sponding to each A,, there is an eigenfunction y,, = X,,(x) satisfying (1) and 





(2) and these functions form a complete orthogonal set in the interval (0, /). 
Any arbitrary function f(#) defined in this interval may be expanded as 


[Quart. Journ. Mech. and Applied Math., Vol. IX, Pt. 3 (1956)] 
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an infinite series of these functions under conditions similar to thos 


governing ordinary Fourier series, viz. 
a 
f(x) = > a, Xx n(®), (3 
n=0 
f 
where a, = | r(x)f(a)X,,(x) dx. (4 


In general, for a given example of (1), the practical determination of the 
eigenvalues and eigenfunctions is not necessarily an easy matter, and 
solutions are usually obtained by standard methods, such as the use of 
power series, or numerical methods involving finite differences. Other 
methods widely used for the derivation of approximate solutions are based 
on variational principles. The object of the present paper is to suggest ; 
numerical analysis based on the result (3). In this application the function 
f(x) is itself the required eigenfunction satisfying the particular equation 
of type (1) which we require to solve, while the functions X, (x) are known 


} 


| 





eigenfunctions of another equation of type (1), each X,, satisfying the | 


conditions (2). The functions X,, will depend upon the type of equation 
(1) to be solved, and it will be shown that (1) may be reduced to an infinite 
set of simultaneous algebraic equations for the coefficients «,, which are well 


> 


conditioned for solution by numerical methods. It will also be shown that | 


the method is basically equivalent to an application of the Rayleigh-Ritz 


variational method, provided the correct variational functions are chosen. \ 


It should therefore be pointed out that the basic principles may not b 
new in special cases, but no general attempt at an exact formulation of the 
problem along the lines suggested here seems to have been made, and thi 
method certainly does not appear to be in use as an accepted numerical 
technique. 

In the present paper we shall consider only the case p(2) = 1. In this 
vase the functions X, are conveniently chosen as the trigonometric fune 
tions which are the solutions of 


X”"+a7X = 0 (a real) (5 


subject to the conditions (2). These include as special cases the ordinary 


Fourier sine and cosine series in the interval (0,1) which arise for particula! 4 


values of h, and h,. Two numerical examples, both derived from physical 
problems, are treated in detail; these have been chosen to give some 
indication of typical numerical computations rather than for any special 
difficulties which they present. 
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It is perhaps relevant to the analysis of this paper to mention here a 

special example of (1), namely, Hill’s equation 

y'+w(xr)y = 0, (6) 


where w(2) is an even periodic function of x given by 


a 
~I 
~— 


x 
w(x) = 0,+2 > 86,,cos 2rx, 
r=1 


where > @,,| converges, and the known constants 6,, may or may not 
include a characteristic parameter. This equation includes Mathieu’s 
equation as a special case. Much literature has been published on the 
solution of (6), of which McLachlan (2) has given a connected account. 
The solution is obtained in general form by assuming a trial solution 
y= or > cyeme, (8) 
pono 
and on substitution in (6), an infinite set of simultaneous equations is 
obtained for the c,,. The admissible values of » are obtained as the latent 
roots, Only two of which lead to distinct solutions, of an infinite deter- 
minantal equation, and the corresponding c,, are obtained numerically 
from the simultaneous equations. The method is quite distinct from 
that proposed here, which, however, includes Hill’s equation as a special 
case: but the important point is that Hill’s determinant and the associated 
simultaneous equations are very similar to those which we shall obtain, 
and much of the discussion which has been given on questions such as 
convergence is equally applicable. 
2. Method of solution 
The eigenfunctions X,, satisfying (5) are given by 

X, = sina,z+p,cosa,2% (n= 0, I, 2,...), (9) 
and these satisfy the conditions (2) provided 
y, == fo. (10) 

x, (Ay —h,) 


fe A 1] 
v2 +h, he ome 


and tana, | 


For specified values of h, and h, the positive roots of (11) define the 
characteristic values of « corresponding to the eigenfunctions X,. These 


functions then satisfy the orthogonality properties 


| X,A,d%&=0 (mn) (12) 


and ; X? da 5 (1 4 Xn fo (hy hs) . hy hs) (13) 
d 2 hi} 2hi(a® +h) 
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With p(x) = 1 in (1), the equation to be solved is 
y” +[q(x)+Ar(x)]y = 0. 
In accordance with (3) we now write 
Yy = = a, Xx, 
n=0 


so that, on using (12), we obtain 


l 
| yX,, dx = 8,4, 
0 
4 
and thus | y”X, dz = —a25,a,, 


0 


on integrating twice by parts. Hence, multiplying (14) by X 


grating with respect to x over the range (0,1) we obtain 


l 


—a2 $a, 4 ( [g(x)+Ar(x) |X, ydx=0 (n= 


0 


It now remains to express the integral in (17) in terms of the coefficients 
a, and, when this has been achieved, (17) will be reduced to a set of simul- 


0), 


i, 


taneous algebraic equations for the a, which we may write 


x 
J ; - m 
> Orm@m =O (n= 0,1 
m=0 


o Mi yeee 


). 


Megeee 


( 
( 
( 


(17 


) 
14 > 
t 
16 | 


» and inte- 4 


a 


(18) 


: ! 
l'o evaluate the terms 5,,,, we substitute for y in the above integral by 


using (15) and note that the functions X,, may be expressed in the form 


X,, = (1+ph)! sin(a, x+B,) | 


where tan fp 


We may now integrate term by term. Thus let us adopt the notation 


l 


au" . Pn J 


Cmin q(x)cos} (x, T X,)& T (Bn»+B, } dx 


0 


and denote also by c,,_,, the value of the integral (20) when the signs of 


x, and 8, are changed. Two similar quantities d 


tities thus defined, it is easily verified that 


Diss = — = 3(1 T pr) + Din) *{Cm, n—Cmnt+A(d 


«ol 2\Sn 
Diss ii 3(1 T PrCn,—n—Can +A(d d 


n,n 


where Cnn = | Q(x) dz, 


m,n 


1 
nn ) f 


m, 


—T 


—s 


and d 


m,-—n 


-d 


a 


mn ) 


n 


j 


(19) 


(20 


denote the 
values of (20) when q(x) is changed to r(x). Then, in terms of the quan- 
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for all values of nm, and d, _,, is similarly defined. The formulae (21) com- 
pletely define the simultaneous system (18), which will have a solution in 
which the coefficients a, are not all zero provided that the infinite deter- 
minant formed by the b,,,,, is zero, that is 


A(A) = ||8am|| = 9. (22) 


i 7m 
The roots of (22) define an infinite number of eigenvalues A,,, provided 
that these roots exist, and this will be the case if it is possible to make 
A(A) convergent. An infinite determinant converges if the infinite product 
of its leading diagonal elements and the infinite sum of its non-diagonal 
elements are absolutely convergent.¢ Thus, dividing each row of A(A) 
by its corresponding leading diagonal element, it will be seen that A(A) con- 


verges prov ided 
>. > & 


m=0 n=0 


(m =~ n) 


nm 


converges, and A has no value which makes a leading diagonal element 
zero. The functions q(x) and r(a) will have to be such as to satisfy these 
conditions. We may then expect that the system of equations (18) will 
be well-conditioned for solution by numerical methods, since ultimately 
the term in a2 must dominate the leading diagonal. Moreover, in some 
problems it is possible that a good first approximation to the eigenvalues, 
particularly the lower ones, may be obtained by setting each leading 
diagonal element successively equal to zero. 

It should be pointed out that the above analysis is equivalent to a 
variational treatment using the Rayleigh—Ritz principle. For the con- 
dition that the equation (14) shall arise is that the expression 

; : 

J | {y"+-q(x)y}y dx +A | r(ax)y? dx 

0 0 
shall have a stationary value. If we substitute in these integrals the value 
for y in terms of the X,, given by (15) we obtain J as a function of 

, Ay yeeey a, genes 

and the condition @.J /éa, O (n 0, 1, 2,...), gives rise to precisely the 
same set of equations defined by (18) and (21). The Rayleigh formula for 
the eigenvalues, which may be used to estimate values of A at any stage 
of a numerical computation, is given by 


} 


l 

” 9 

m\ | fy” +q(a)yhy del r(a)y* dx, 
0 0 

+ A brief account of infinite determinants is given by Whittaker and Watson (3); the 

same authors discuss Hill’s determinant and the associated simultaneous equations. See 


also the remarks of Jeffreys (4) on infinite determinants, and the footnote to p. 382 of the 


present paper 
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so that, in terms of our present notation, we obtain the formula 


x x x 
S 2 72 1 y 2\3 2 \h » 
> 0, Xn an—%s > e (J T Pr)? (1+Pm)? (Cm, n —Cmn)OBp Am 
j= *% 0 __ n=0 m=0_ paw ees _ (23 
x x “v0 
1 2\4 2 \} 
2 > > (1+ ae % Pin)* (din, —-n— Finn On am 
n=0 m=0 


where, in the double summations, the product of each pair of non-diagonal 
terms a, a,, occurs twice with equal coefficients, due to symmetry about the 
leading diagonal. In many practical problems the formula (23) and the 
previous analysis will simplify, on account of special forms of g(x) and r(z’, 
For example, in some classes of problems of elastic stability g(x) = 0, 
while in many problems in which q(x) 4 0 it is common for r(x) to be 
constant, a case in point being Mathieu’s equation. 


3. Special forms of the boundary conditions 

For sufficiently large values of n and for all finite, non-zero, values of 
h, and hg, we find from (11) that «,, ~ nz/I so that ultimately the coefficients 
Cm» and c,,_, become the members c,,,,, and c,,,_,,, of the singly infinite 
set of coefficients of the Fourier cosine series 


x 


c 7 on Nrxe 
q(x) = 7+ ; 2 t. cos —- ; (24) 


n=1 

Similar remarks apply to the coefficients d,,,, and d,, _,,.. To the accuracy 
required for most computational purposes, this result applies for even 
moderately large values of n and this facilitates computation of the 5,,,.. 
For certain forms of the boundary conditions, however, which include 
those most commonly occurring in practice, this result holds for all values 
of n and the whole analysis may be performed in terms of ordinary Fourier 
series in the interval (0,1). We now briefly consider these cases. 

(i) y(0) = y(l) = 0. In this case we put h, = hy = & so that p, = 0, 
5, 41 from (13), and 8, = 0 from (19). The roots of (11) become the 
positive roots of 


sin «,, / 0, (25) 

nt : a 

that is X » ; (ns = 0, I, 2.,...), (26) 
, , - NTL i 

so that A, = 0, XxX, = sin j (n = |, 2, 3....) (27) 


and y is represented by a Fourier sine series in (0,/). If we put 8, = 
in (20), the coefficients defined by this equation all become members of 
those defined by the series (24), and corresponding to (21) we obtain the 


results Ls 
b, m Iman 21Cim—n Cmin A(d m—n| dn n ) j | 
‘ » 
1 | n>?) . (25) 
Inn 2 \° Con A(dy dy,,)— > l { 
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(ii) y/(0) = y'(l) = 0. Let us first multiply each equation of the system 
18) by A} and write a5 5;,, = Aj5,. We may now put h, = h, = 0, so that 
3, = 4, while from (19) we have B 


11) are again given by (25) and (26) and the eigenfunctions X,, are now 


= 4m since p, = 0. The roots of 


a 


X, cos—_— (s = 0, I, 2....) (29) 

l 
so that y is given by a Fourier cosine series in (0,1). On putting 8, = 47 
in (20) we obtain c,,, ,, Cn+n And C,, , = Cm»), With a similar result 
for d,,,, and d,, ,,, so that the 6,,,, are derived from (28) with appropriate 


changes of sign. 


(iii) y(0) = 0, y’(l) = 0. Here we put h, = 00, h, = 0 so that as in (i) 


we have p 0 and £, 0. The roots of (11) are now the roots of 
cos x, l c= @, 
, ni a" : 
that is Y 7 (n ES. S...a (30) 
, . NX at. 28 " 
ind X, sin (n = 1, 3, 5....). (31) 
2] 
Thus y is represented by the Fourier sine series in (0, 2/) which is sym- 
metrical about x = /. Again in this case we may express the b,,,,, entirely 


in terms of the coefficients of the series (24) and the corresponding series 
for r(x). Thus, subject to the restrictions on n in (30) and (31), we obtain 


Lig » im 
b, m Moss 2 uf 5 m—n c Mm+n) T A(d, mi diimsn))} | 
n2z2 7” (32) 
Ban = 5\¢o—Cu + Ady—d,) — "| | 
2 | 4] | 


In all cases of the previous analysis it should be possible, for many 
analytical forms of g(a) and r(x), to derive simple formulae for the 5,,,,, in 


nm 
terms of a, and «,,. If this is not the case, or these functions are specified 
numerically, the process is still feasible by using numerical integration and 
in this event the special integration formulae of Filon (5) may prove useful. 


These formulae are, in effect, modified Simpson formulae for trigono- 
b 

metric integrals of the form | (x)sin(ka+e) dx which, by proper choice 
a 


of k and e, can be used to evaluate (20). They are not given in the present 
paper since they are most efficiently used in conjunction with the table 
of appropriate factors given in Filon’s paper. 

The remainder of this paper is devoted to the solution of two numerical 
examples. Neither presents a specially difficult application of the method, 
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but both are derived from physical problems, and our main object is to 
indicate that in problems of this kind the numerical technique used jp 
their solution can be quite rapid. 


4. A problem in quantum mechanics 
We shall first consider the numerical solution of a well-known problem 
in the field of quantum mechanics, namely, that of determining the 


vibrational wave functions of the one-dimensional harmonic oscillator. 
In this case the Schrédinger equation for the wave function ws is (see (6)) 


» 
a ee 


4+—~(E—V)p = 0, 
h2 
where the potential function V = ax?. If we put 2ma/h? = o2, x = x 14 
this equation becomes ‘ we - 
{Ue pb” +(A—L2)b = 0, (33) 
E(2m\4 : : ; 
where A = j “. If the oscillator is unbounded the boundary con- 
r\a i 


ditions are that %(co) = y(—0oo) = 0, and the solutions of (33) are 


. = | 
ian ( eH (0) (34) 


(2" n!)* \a 


n 


provided that A, = 2n+1 (n = 0,1, 2....), (35) 


where /,,(£) is the Hermite polynomial of order n, so that the wave functions 
(34) are even or odd with n. 
The eigenvalues (35) correspond to the energy states 


’ 


EB, = (n+4)ho, 


where w {= (2a/m)'} is the angular frequency of vibration of the oscillator. 
In the case of the bounded, or enclosed, oscillator, of some importance in 
diamagnetism, the boundary conditions are that % shall vanish at ¢ = +4 
where }/a~* is the amplitude of the oscillator. This problem has been con- 
sidered by Auluck and Kothari (7) inter alia; for given finite values of | 
they determine the energy levels ,, by investigating numerically the zeros 
of the appropriate Hermite functions, the results being given graphically. 
It is also shown that exact expressions for the energy levels may be ob- 
tained in the two limiting cases of / sufficiently small or sufficiently large. 
The investigation for small / is of particular interest in the present paper: 
in this case they obtain the expression 

7h? pax 


gE. = 


Pp 


= 1, 2, S.... 36) 
2ml? (P ) ( 


which is derived from the limiting form of the Hermite functions. A second 
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approximation is obtained by substituting in (33), after first changing the 


ange of integration to (0,/), the solution 
rang 


b= S Asin _ (a == 1, % &...., (37) 
m=1 


which is equivalent to our form of solution. No attempt is made, however, 
to formulate the simultaneous equations, but approximate analytical 
formulae are derived for the eigenvalues and eigenfunctions.t 

The problem presents an application of the formulae of § 3(i). If we 
change the origin by the substitution = ¢’—4I, we obtain 

q(C') (¢’—3l)? and r(f’)=1. 

The odd coefficients c,, are found to be zero so that the non-diagonal terms 
b are zero unless m and n are both even or both odd. Subject to this 
restriction we find 


3 } 2.2 
h 4mnl . b ; If 1p3(1. 6 mL r*| 
2 | ? n J 


a*(m*—n?)? ‘sis 2772 


and the simultaneous system of equations (18) breaks up into two inde- 
pendent sets, one each for the odd and even coefficients a,. Similarly, 
the infinite determinant A(A), which is seen to satisfy the conditions for 
convergence, factorizes into two determinants, one each for the even and 
odd states respectively (note that the pth eigenfunction in our notation 
corresponds to the (p—1)th characteristic state as defined in the notation 
of (34)). Approximate values of A may be obtained by expanding the 
infinite determinants about the top left-hand corner and taking into 
account sufficient columns and rows, although this method soon becomes 
laborious if more than a few columns and rows are required. Some results 
are given for the even states in the three cases 1 = }2, 7, and 27 in Table 1, 
where for a given value of / the four approximations correspond to deter- 
minants of the first, second, third, and fourth order respectively. As / is 
increased the eigenvalues tend to (35). 

Probably the best method of computing the eigenfunctions and eigen- 
values, however, is to solve the simultaneous equations by standard 
relaxation or iterative methods. For example, the relaxation method 
can be used in conjunction with the use of Rayleigh’s principle or the 
method of intensification.t In the present example, and the one of the 
next section, the following iterative method was found to be very rapid. 
Supposing that we wish to compute the eigenfunction %,, and its associated 
eigenvalue A,,, we first put a, = 1 in each equation of the system (18) and 

+ 


[t should be noted that equation (19) of the paper cited is only an approximate result. 
{ An account of the practical use of these methods has recently been given by Allen (8). 
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write the remaining non-diagonal terms of each equation equal to § 
n 


(n = 1, 2, 3,..., p,..), so that these equations become 


—(b,,,+6 
a, = bh -— (n a Pp), (38) 
[2 6 pn? = 28 
eee 3b. Oe Ae 3$ 
p l (! ani ]2 l \ se) 


TABLE | 





R | Approximation A; As As Ay 
1st 4°0806 
1 2nd 4°0803 | 36°192 
5: 3rd 4°0803 | 36°192 | 100°20 
qth 4°0803 | 36°192 100°20 196°20 
1°3225 
*- 1*3055 9°7835 
1*3057 9°7697 25°816 
1°3057 9°7695 25°806 49°822 
1°5399 
—_ 1°0167 5°8408 
10004 5°0587 10°258 
1°0003 5°0464 9°6233 16°146 





Initially we assume all the 5, zero and then, using A,, calculated from 
(39), we find estimates of the a, from (38). From these the 8, may be 
calculated and fresh estimates of \,, and the a,, determined and the process 
repeated until complete convergence is achieved. At any stage the Rayleigh 
estimate of A,, given by (23) will be more accurate than (39), but the latter 
was found adequate in the present calculations. The results of some 
typical calculations, correct to an overall accuracy of five significant 
figures, are given in Table 2. In the cases considered the first estimate of 
A,, from (39) was quite good. The formula (36) can of course be obtained 
from (39) for, when / is sufficiently small, we find A,, ~ p?x2/l?, which is 
equivalent to (36); a more accurate formula can be obtained by using the 
full formula (39) with 6, = 0. 
TABLE 2 





R Aor R 7 | R 27 
A; 4°0803 | A, 1°3057 | As = 9°7694 | Ay = 170003 
1 Im am am Im 
I 1*0000 10000 | O'0441 1*0000 
3 0°0029 0°0442 I-000C 0° 3687 
5 "0002 00020 | 070294 | 0°0495 
7 0°0004 0°0023 | —0*0027 
9 } 0*0001 0°0005 | 
Ir | | 00002 
| 
| 
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5, A problem in heat transfer 

As a second numerical example we shall consider the problem of the 
determination of the temperature in a viscous incompressible fluid moving 
with Poiseuille velocity distribution between infinite parallel plates; the 
(uid enters the plates at temperature 0, and the plates are capable of 
finite heat transmission into a medium at zero temperature. Taking x and 
y as longitudinal and transverse coordinates with respect to the plates, it 
may be shown (see, for example, (9)) that, subject to certain approxima- 
tions, #(2, y) is the solution of the equation 


0 — 3Pé(1—4y?)‘ (40) 
cy* 4 Ox 
with a) 6, x= 0, _ 4 = y “ 4 
c0 6 (41) 
hd aty i, - =—h? aty=}4,27r>0 
cy - cy i 


where Pé and h are thermal constants, and the plates have the same thermal 
properties. The problem may be solved by substituting in (40) the solution 


6=YSA 
n=1 


e -2Azr/8PEY (y), (42) 


where the functions Y, (y) are eigenfunctions of 
Y”"+A(1—4y2)¥ = 0 (43) 
with Y’(—4) hY(—}), Y'(4) = —AY (4). (44) 


The problem has in fact been solved, for several values of h, by van der 
Does de Bye and Schenk (9) who obtain a power series solution of (43), viz. 


x 
_ 2r 
Y = Yay - 

n=0 


where 2(n+1)(2n+-1)a, ,,+A(a, —4a,_,) = 0. (45) 


Two points arise in connexion with (45): (i) the terms a,, increase as A" and 
a large number of terms of the series may be required when A is large 
(cf. A, = 293), and (ii) each eigenvalue is obtained by tabulating Y and 
its derivative at « = } and determining the value of A for which (44) is 
satisfied, a process which may not be very sensitive. 
In the present mode of solution we change the origin by the substitution 
; 


y = »—4 and then, in terms of previous notation, we have q(n) = 9, 


/ 


r(7) 4n(1—7n), 1 |,andh, = —h, = h. The roots of (11) now become 
the roots of 7 
x, tanta, = h. (46) 


The first six roots of (46) have been tabulated by Carslaw and Jaeger (10) 


for a wide range of variation of h. The terms d d 


m,n? m,—n 


are obtained by 
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integration and, after a little simplification, we obtain the results 
= 2 2 2}2)/ 2 
a = —8Aa,», X,, ( Xm TT me 4h-+-2h? )/ h? (a2, —t ~ 4 


‘| x2 2 x2 —h? l x2 ] 
b _— A l Ci) a oe x bs Be nm) 
- 3 = rh AF | 3/5 ( | 3) rh 


Computation of the first three eigenfunctions has been carried out in the 
case h = 2 by use of the method described in the previous section. Results 
have also been given in this case by van der Does de Bye and Schenk, 
and a comparison of the first three eigenvalues is given in the following 
table; the quantity tabulated is B,, = (A,,)!, which they tabulate. 


.T 





By Bo Bs 
van » dee Does de Bye and Schenk , . | 20000 9°3124 
Present results 


24 
24 


I 
I 


I 
2*0000 | 9° 3123 I 


| 


TABLE 3 








A, = 4°0000 | A, = 86°719 | As = 293°23 
" Im am Im 
I 1*0000 o'8oll I*2119 
2 O°O105 1*0000 1°1864 
3 0°0003 o1116 1*0000 
4 0°00 30 0°3123 
5 0°0009 + 0°0057 
6 0°0003 0°0027 
7 0°0001 o*0010 
5 0°0001 0°0005 
9 0*0002 
10 0°000I 
11 


o°ooo! 


Our results for the eigenfunctions are given in Table 3. Having determined 
sufficient of the eigenfunctions, the coefficients A,, in (42) may be obtained, 
in the usual way, by making use of the first of (41) and the orthogonality 
properties of the functions Y,. Explicit formulae may be derived in terms 
of the a, following the manner described in (1), in which the case h = 0 
of the present problem has been solved. In general the evaluation of these 
coefficients, and also of relevant physical data, is greatly facilitated by the 
fact that the orthogonal functions Y,, are expressed in terms of the functions 
X,,, which are themselves orthogonal in the interval (0,1), so that simple 
formulae in terms of the a, can usually be obtained. For example, it can 
be shown that the constants A, are given by the formula 
A,, = 20, : oy a, > Y Sno @e, (n = 1,2,3...), 


m 
m 
+ We must point out Ps in this actin the terms by» do not satisfy the conditions 
of convergence of A(A) stated on p. 375. 
ever, 


This apparent difficulty may be removed, how- 
by writing «3,4, = aj, in (18) and forming the determinant A’(A) for non-vanishing 
Qm, Which is convergent. 
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THE 
where each value of n corresponds to the appropriate set of a,, given in 

Table 3. 

* Inconclusion, it should be re-emphasized that the examples treated in this 
and the previous section do not present difficult applications of the method, 
though such could doubtless be constructed. They have been chosen 
because they are typical for many analytical forms of g(x) and r(x). For 
example, for many of the elementary functions the coefficients c,, and d,, in 
the series (24) and its counterpart will decrease like 1/n?, and in such cases 


the convergence of the numerical process should be rapid since the large 


this is not the case, however, the computations are likely to be more 


f 
elemenis 6,, will be near the leading diagonal. In any problem in which 
difficult. 
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